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ABSTRACT 


High  heat  flux  removal  from  devices  such  as  Insulated-Gate  Bipolar  Transistor  (1GBT) 
and  Monolithic  Microwave  Integrated  Circuit  (MM1C)  will  be  important  for  Navy  ships 
in  future  applications.  Micro  encapsulated  phase  change  material  (MEPCM)  slurry  was 
used  as  the  heat  transfer  fluid  inside  a  microchannel  instead  single  phase  fluid  in  present 
study.  Presence  of  phase  change  material  increases  the  effective  heat  capacity  of  the 
fluid.  To  enable  the  application  of  using  MEPCM  slurry  in  microchannel  cooling,  three 
research  tasks  were  identified;  namely,  clogging  experiments  (Task  1),  heat  transfer 
experiments  (Task  2)  and  numerical  modeling  (Task  3).  While  prior  reports  have 
described  the  efforts  in  Tasks  1  and  2  and  the  corresponding  results,  the  focus  of  this 
report  is  on  Task  3  involving  numerical  modeling  of  EPCM  slurry  in  microchannels. 

Particle  migration  inside  microchannels  of  different  aspect  ratios  was  simulated  using  the 
lattice  Boltzmann  method  (LBM).  Channels  of  aspect  ratios  from  1  to  4  have  been  used 
for  modeling  the  geometry.  The  particle  size  to  the  smallest  dimension  of  the  channel 
ratio  is  taken  to  be  1/10  and  the  particle  concentration  is  less  than  0.5%  for  all  cases 
simulated.  Simulations  show  that  these  particles  migrate  to  preferred  locations  in  the 
channel,  known  as  equilibrium  locations,  and  these  depend  on  Re  and  aspect  ratio  of  the 
channel.  For  aspect  ratio  2  and  higher,  particles  form  an  inner  ring  that  moves  away  from 
the  wall  as  Re  is  increased  further.  An  outer  ring  is  formed,  which  moves  closer  to  the 
walls  as  Re  is  increased.  The  location  of  these  inner  and  outer  rings  is  strongly  dependent 
on  the  Reynolds  number  of  the  flow. 

For  modeling  higher  concentrations,  a  shear  induced  migration  model  based  on 
macroscopic  constitutive  equations  was  employed  and  was  used  to  simulate  the 
nonhomogeneous  particle  distribution  inside  the  microchannel.  The  particle  distribution 
was  modeled  using  the  diffusive  flux  model  and  the  results  were  used  to  solve  the 
temperature  profile.  In  order  to  investigate  the  effect  of  particle  migration,  the  thermal 
results  obtained  with  and  without  including  the  particle  migration  were  compared.  It  was 
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found  that  for  the  parameters  considered,  the  difference  in  slurry  performance  prediction 
assuming  nonhomogeneous  and  homogeneous  distribution  is  not  significant. 

For  parametric  study,  the  three  dimensional  flow  inside  the  microchannels  was  solved 
including  the  microchannel  wall/fin  effects  and  developing  flow  in  manifold 
microchannels.  A  constant  inlet  velocity  and  temperature  were  assumed  along  with 
constant  heat  flux  condition  at  the  base  of  the  microchannel.  Thermal  performance  of 
water  and  poly-alpha-olefin  (PAO)  based  slurries  were  analyzed  in  microchannels  of 
width  101  pm  and  25  pm.  Pressure  drop  of  both  slurries  is  higher  compared  to  the  pure 
fluid  as  expected.  Heat  transfer  coefficient  of  slurry  is  lower  compared  to  water  because 
of  lower  thermal  conductivity  in  microchannels  of  w  idth  101  pm.  In  case  of  25  pm  w  ide 
channels,  water  based  slurry  performed  better  compared  to  water  due  to  faster 
development  of  temperature  profile.  It  was  found  that  PAO  slurry  performs  better  in  both 
channels. 
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PROGRESS  STATEMENT  SUMMARY 


The  performance  of  encapsulated  phase  change  material  (EPCM)  slurry  flow  in 
microchannels  was  investigated  using  the  effective  specific  heat  capacity  method.  Lattice 
Boltzmann  method  was  used  to  simulate  the  particle  paths  when  the  duct  shape  has 
different  aspect  ratios.  For  higher  concentrations,  a  shear  induced  migration  model  was 
used  to  simulate  the  nonhomogeneous  particle  distribution.  Results  of  the  model  were 
used  to  solve  the  temperature  distribution  inside  the  channels.  Parametric  study  was 
continued  with  water  and  PAO  as  base  fluids  in  microchannels  of  two  different  widths, 
101  pm  and  25  pm.  Parametric  study  was  done  by  varying  parameters  such  as  particle 
concentration,  inlet  temperature  of  the  fluid,  melting  range  of  PCM,  base  heat  flux  and 
base  fluid. 
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a 
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Bip 
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Aspect  ratio 
Area  of  base,  m* 

Particle  radius,  m 

Width  of  the  channel,  diameter  of  the  tube,  m 

Biot  number 

Biot  number  of  particle 

Volume  concentration  of  MEPCM  particles  in  slurry 

Mass  concentration  of  MEPCM  particles  in  slurry 

Specific  heat  of  bulk  fluid,  J/kg.K 

Specific  heat  of  fluid,  J/kg.K 

Specific  heat  of  MEPCM  particle,  J/kg.K 

Specific  heat  of  PCM,  J/kg.K 

Hydraulic  diameter,  m 

Particle  diameter,  m 

Magnitude  of  shear  rate.  I/s 

Phase  space  vector  along  link  "i",  m/s 

Enhancement  factor 

Particle  distribution  function 

velocity  gradient,  s'1 

effective  velocity  gradient,  s'1 

Height  of  the  channel.  Shortest  dimension  of  rectangular  channel,  m 

Height  of  manifold  in  simulation  domain,  m 

Heat  transfer  coefficient,  W/m‘.K 

Heat  transfer  coefficient  ratio,  =  hsiurry/hbase  fluid 

Latent  heat  of  fusion,  J/kg 

Thermal  conductivity,  W/m.K 

Thermal  conductivity  of  bulk  fluid,  W/m.K 

Effective  thermal  conductivity  of  bulk  fluid,  W/m.K 

Thermal  conductivity  of  fluid,  W/m.K 
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m 
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Pe 
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q 

qw 
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Re 
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Rep 

rp 
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r 

rP 
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T, 

T2 
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Thermal  conductivity  of  MEPCM  particle,  W/m.K 
Thermal  conductivity  of  PCM,  W/m.K 
Length  of  the  channel,  m 

Length  of  the  channel  for  fully  developed  particle  distribution  profile,  m 

Mass  flow  rate  inside  the  microchannel,  kg/s 

Nusselt  number 

Unit  normal 

Peclet  number 
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Prandtl  number 
Pressure,  psi 
Pressure  at  the  inlet,  psi 
Pressure  at  the  outlet,  psi 
Phenomenological  constants 
Heat  supplied,  W 
Heat  flux,  W/cm" 

Constant  wall  heat  flux,  W/cm 
Radius  (m) 

Rey  nolds  number 

Reynolds  number  for  a  circular  channel 

Particle  Reynolds  number  =  Remaxfa/H)2 

Radius  of  MEPCM  particle,  m 

Space  coordinate 

Radial  coordinate,  m 

Solid  liquid  interface  radius,  m 

Temperature,  K 

Lower  melting  temperature.  K 

Higher  melting  temperature,  K 

Temperature  at  heat  sink  inlet.  K,  Temperature  at  microchannel  inlet.  K 
Melting  range  =  T1-T2,  K 
Melting  temperature,  K 
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Tout  Temperature  at  heat  sink  outlet,  K,  Temperature  at  microchannel  outlet,  K 

Tw  Wall  temperature,  K 

Twmax  Maximum  wall  temperature,  K 

t  Time,  s,  Time  coordinate 

tbase  Base  of  microchannel,  m 

tuh  Half  of  microchannel  w  idth,  m 

tres  Residence  of  MEPCM  particle  inside  the  channel,  s 

4  Stokes  time,  s 

tw  Half  of  MicroChannel  wall  thickness,  m 

U  Velocity  distribution,  m/s 

UaYg  Average  velocity  of  the  flow,  m/s 

Ueff  Effective  velocity  on  the  particle,  m/s 

Uo  Maximum  velocity  of  the  flow,  m/s 

u,  u  Velocity  vector,  m/s 

u  Velocity  in  x  direction,  m/s 

Vhs  Volumetric  flow  rate  at  the  heat  sink  inlet,  m3/s 

Vr  Radial  v  elocity  of  the  particle,  m/s 

Vz  Axial  velocity  of  the  particle,  m/s 

v  Velocity  vector 

v  Velocity  in  y  direction,  m/s 

W  Half  channel  width,  m 

w  Velocity  in  z  direction,  m/s 

x,y,z  Spatial  coordinates 

x  Distance  from  channel  symmetry,  m 

x  Space  vector 

Y  Non-dimensional  coordinate  along  the  *y’  axis 

Z  Non-dimensional  coordinate  along  the  'z’  axis 

AP  Pressure  drop,  psi 

APr  Pressure  drop  ratio,  APS|Urry/APbasetiuid 

ATbuik  Bulk  temperature  rise,  T0ut-Tjn,  K 

At,  dt  Time  step,  s 
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Sx 


Axial  relative  velocity  of  the  particle,  m/s 
Grid  spacing,  m 
p  Pressure,  Pa 

Pe  Peclet  number 

Greek  symbols 


a 

Thermal  diffusivity,  m2/s 

ap 

Thermal  diffusivity  of  MEPCM  particle,  m2/s 

ar 

Thermal  diffusivity  of  fluid,  m  / s 

Y 

Shear  rate,  1/s 

K 

Velocity  gradient,  s'1 

0 

Non-dimensional  temperature  =  (T  -  Tm)/  (qw  Rd/kbuik) 

p>n 

Dynamic  viscosity  (Pa.s) 

Pb 

Dynamic  viscosity  of  bulk  fluid.  Pa.s 

pt 
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pP 
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V 
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p 
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pb 
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pf 

Density  of  fluid,  kg/m‘ 

PP 

Density  of  MEPCM  particle,  kg/m' 

r 

Time  constant,  s 

<P 
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b 
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d 
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eff 

Effective 
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eq 
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f 
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/ 
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/ 
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m 
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max 

Maximum 

P 
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SECTION  1:  TECHNICAL  OBJECTIVES 


We  have  identified  three  key  fundamental  issues  in  the  application  of  micro  encapsulated 
PCM  slurries  in  microchannel  heat  sinks.  Three  tasks  are  proposed  which  will  improve 
our  understanding  of  the  issues,  resulting  in  solutions  to  potential  problems.  There  have 
not  been  any  experiments  reported  so  far  that  use  PCM  slurries  in  channels  of  less  than 
3.14  mm  size.  Hence  it  forms  an  interesting  and  important  aspect  to  perform  experiments 
on  microchannel  heat  sinks  with  channel  sizes  between  100  pm  and  500  pm  and  filled 
with  MEPCM  slurry  with  different  particle  sizes  and  concentrations.  Also,  the  validity  of 
the  existing  modeling  work  is  questionable  and  further  theoretical  investigation  is  needed. 
Two  different  experiments,  one  for  observing  the  clogging  and  the  other  for  finding  the 
pressure  drop  and  thermal  performance  of  slurry  will  be  performed.  The  numerical 
modeling  will  include  the  modeling  of  particle  distribution  in  microchannels  and  thermal 
performance  based  on  the  non  homogeneous  distribution. 
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SECTION  2:  TECHNICAL  APPROACH 


For  observing  clogging  inside  the  microchannel  a  microfluidic  channel  test  section 
fabricated  at  UCF  were  employed.  These  channels  are  100  pm  in  height,  500  pm  in  width 
and  6  mm  in  length.  Water  will  be  used  as  base  fluid  and  MEPCM  particles  with  n- 
eicosane  as  phase  change  material  will  be  used  for  experiments.  This  experiment  will  be 
done  with  particle  diameters  less  than  30  pm  and  for  different  Reynolds  numbers. 


For  the  heat  transfer  and  pressure  drop  experiments,  a  custom  fabricated  SA-2  cooler 
from  Micro  Cooling  Concepts,  Inc.®  will  be  used.  These  channels  are  100  pm  in  width, 
500  pm  in  height  and  I  mm  in  length.  In  these  experiments,  heat  transfer  and  pressure 
drop  across  a  SA-2  cooler  using  MEPCM  slurry  as  the  working  fluid  will  be  analyzed 
under  different  concentrations  and  Reynolds  number  and  the  results  will  be  compared 
with  those  of  using  water  in  the  same  test  setup. 


Since  MEPCM  slurry  flow  in  microchannels  is  a  kind  of  suspension  flow  of  particles  in 
fluid  under  motion,  the  particle-particle,  particle-wall  and  particle-fluid  interactions  are 
important.  Lattice  Boltzmann  approach  is  a  realistic  method  to  simulate  suspension  flows 
as  it  accounts  for  all  the  forces  is  not  as  computationally  intense  as  molecular  level 
dynamic  simulations.  The  lattice  Boltzmann  method  (LBM)  is  based  on  a  microscopic 
model  but  is  used  to  model  macroscopic  fluid  properties.  This  method  is  used  for 
simulating  the  particle  distribution  inside  the  microchannel.  A  macroscopic  model  that 
simulates  the  particle  distribution  along  the  channel  based  on  the  diffusive  flux  model 
will  be  used  for  simulating  the  particle  distribution  inside  the  microchannels.  These 
results  will  be  used  for  analyzing  the  thermal  performance  of  slurry.  For  this  purpose, 
using  the  simulation  results  shear-induced  hy  drodynamic  model  for  the  MEPCM  slurry  in 
microchannels,  a  volume  averaged  macroscopic  model  on  heat  transfer  will  be  developed 
based  on  continuum  energy  transport  equation.  The  effective  thermal  conductivity  will  be 
calculated  based  on  the  concentration  of  the  MEPCM  particles.  Heat  absorption  delay 
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will  be  considered  by  solving  the  heat  transfer  process  in  the  particles.  The  solution  will 
consider  the  melting  temperature  range  and  the  possible  effect  of  the  sensible  heat  of 
PCM  on  heat  transfer.  Heat  transfer  coefficient  will  be  calculated  based  on  the  liquid  and 
particle  velocities  obtained  from  the  hydrodynamic  simulation.  Heat  absorption  will  be 
considered  as  a  heat  sink  term  in  energy  equation  for  liquid.  The  analysis  will  be 
compared  with  the  experimental  results. 
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SECTION  3:  ANNUAL  PROGRESS 


3.1  Introduction 

Tuckerman  and  Pease  in  1981  demonstrated  a  microchannel  heat  sink  that  removes  790 
W/cirf  with  71  °C  temperature  increase  at  600mL/min  flow  rate  [1].  The  heat  sink  was 
numerous  small  channels  and  fins  arranged  in  parallel  directly  fabricated  in  a  silicon 
substrate  and  had  direct  circulation  of  water.  It  is  so  compact  that  heat  is  efficiently 
carried  from  the  substrate  into  the  coolant  because  of  its  inherently  small  passageways 
and  a  very  large  surface-to-volume  ratio.  It  has  already  been  proved  that  the  use  of 
microchannel  heat  sinks  has  numerous  advantages  compared  to  macroscale  flow 
channels.  Tuckerman  and  Pease  predicted  that  single-phase  forced  convective  cooling  in 
microchannels  should  be  feasible  for  circuit  power  densities  of  more  than  1000  W/cm  . 
However,  the  heat  sink  had  quite  a  large  pressure  drop  of  200kPa  with  plain 
microchannels  and  380kPa  with  pin  fin  enhanced  microchannels.  For  single  phase 
cooling,  the  coolant  temperature  will  increase  in  the  flow  direction  as  it  acquires  heat, 
which  leads  to  non-uniform  temperature  distribution  on  the  chips.  In  order  to  keep  the 
temperature  of  devices  such  as  semiconductors  and  lasers  to  be  cooled  within  a  few 
degrees  Celsius,  a  very  large  mass  flow  rate  is  needed,  resulting  in  high  pressure  drop 
across  the  microchannels,  even  prohibitively  high  sometimes,  which  need  a  large  pump  in 
the  system  and  consume  a  large  amount  of  pumping  power.  Boiling  heat  transfer  inside 
microchannels  has  been  investigated  to  utilize  latent  heat  of  vaporization  [2],  However, 
the  increased  pressure  drop  and  associated  pressure  fluctuation,  and  wall  temperature 
fluctuation  in  microchannels  hinder  the  application  of  convective  boiling  in 
microchannels  in  electronic  cooling.  Possible  dry-out  at  relatively  low  heat  flux 
compared  with  its  single-phase  flow  counterpart  prevents  two-phase  flow  in 
microchannels  to  be  used  in  cooling  of  high  heat  flux  electronics  [3], 


In  applications  where  high  temperature  variation  and  large  pressure  drop  are 
unacceptable,  manifold  microchannel  heat  sink  (MMC)  [4]  is  shown  to  be  an  effective 
way  in  reducing  the  pressure  drop  and  significantly  reducing  the  temperature  variation. 
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This  is  possible  because  of  the  presence  of  multiple  short  channels  in  the  MMC  heat 
sinks  distributed  along  the  length  of  the  dev  ice  that  needs  to  be  cooled.  The  MMC  heat 
sink  (figure  1 )  features  many  inlet  and  outlet  channels  alternating  at  a  periodic  distance 
along  the  length  of  the  microchannel  [5].  The  flow  enters  the  microchannels  from  the 
manifold  inlet  channel,  splits  and  flows  through  the  microchannels,  then  exits  to  the 
manifold  outlet  channel.  This  pattern  is  repeated  along  the  length  of  the  chip.  These 
channels  are  shown  to  have  small  pressure  drop  compared  to  traditional  microchannels 
with  little  effect  on  thermal  resistance. 


Fin  Base 


Figure  1:  MMC  heat  sink  geometry 


The  intention  of  this  work  is  to  cool  the  high-power  devices  that  produce  high  heat  fluxes 
and  require  continuous  heat  removal  under  near  isothermal  condition.  This  project 
proposes  to  make  use  of  the  latent  heat  of  a  microencapsulated  solid-liquid  phase-change 
material  (MEPCM)  as  working  fluid  to  increase  thermal  capacity  of  the  coolant.  PCM- 
filled  particles  will  be  well  mixed  in  the  liquid.  The  use  of  the  two-component  fluid 
would  provide  additional  thermal  capacity  from  the  latent  heat  associated  with  the  solid- 
liquid  phase  change. 
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The  PCM  in  the  microcapsules  can  be  selected  to  melt  and  freeze  at  the  desired 
temperature.  By  microencapsulating  the  PCM,  the  core  material  remains  separated  from 
the  carrier  fluid,  thus  preventing  its  agglomeration  or  deposition  within  a  system.  In  1983, 
Colvin  et  al.  [6]  at  Triangle  Research  and  Development  Corporation  (TRDC)  in  Research 
Triangle  Park,  North  Carolina,  began  a  long  series  of  pioneering  investigations  that  have 
resulted  in  microPCMs  of  high  strength  enough  to  expose  to  high  pressure  and  external 
media.  For  the  first  time,  TRDC  investigators  demonstrated  that  paraffinic  PCMs  could 
be  successfully  microencapsulated,  suspended  in  water  and  pumped  in  a  closed-loop 
system  [7],  Significant  improvements  were  made  by  TRDC  regarding  the  longevity  of  the 
microPCM  coolants,  which  outlasted  the  bearings  in  the  circulating  pumps  through  over 
1 00.000  thermal  cycles  [6]. 


In  MEPCM  slurries,  the  materials  used  for  the  core  and  the  shell  are  selected  compatible 
and  the  shell  material  is  selected  compatible  with  the  carrier  fluid.  Particle  diameters 
typically  range  from  6  pm  to  50  pm.  The  shell  can  be  less  than  one  micron  thick.  Upon 
heat  absorption,  the  core  material  in  the  particles  melt,  changing  from  solid  to  liquid  and 
the  resulting  volumetric  expansion  of  about  12  to  15%  is  accommodated  by  particle’s 
flexible  wall.  Though  the  density  of  the  PCM  is  usually  not  identical  to  that  of  the  carrier 
fluid,  the  thickness  of  the  shell  wall  can  be  adjusted  so  that  nearly  neutrally  buoyant 
suspension  can  be  obtained  with  relative  homogeneity.  Particle  concentrations  (loadings) 
can  range  up  to  35%  solids  (weight/weight)  depending  on  the  nature  of  the  carrier  fluid. 
In  addition,  the  degree  of  uniformity  may  be  increased  by  adding  a  suitable  dispersing 
agent  [6].  Figure  2  shows  the  cross  section  of  a  single  microPCM  particle,  while  figure  3 
illustrates  a  scanning  electron  microscope  image  of  typical  microencapsulated  PCM 
particles. 
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Outer  Shell 


Figure  2:  Cross  section  schematic  of  a 
single  MEPCM 


Figure  3:  SEM  image  of 
microencapsulated  PCM  Particles 


A  sample  of  PCM  slurry  was  made  at  University  of  Central  Florida  with  water  as  base 
fluid  and  octadecane  as  PCM.  The  particle  size  is  around  0.1  pm  and  the  shell  thickness 
is  around  5nm.  It  is  prepared  by  emulsification  process  where  a  surfactant  (shell  material) 
is  first  dissolved  in  the  base  fluid  and  is  then  heated  to  a  temperature  greater  than  the 
melting  point  of  PCM  while  stirring  continuously.  Later  the  PCM  is  added  to  the  mixture 
to  carry  out  the  emulsification.  The  formation  of  the  encapsulated  PCM  was  confirmed 
using  the  Tyndall’s  effect.  This  involves  the  scattering  of  a  laser  beam  as  it  passes 
through  the  emulsion  due  to  the  presence  of  encapsulated  particles  (figure  4).  Pure  water 
itself  cannot  scatter  the  laser  beam.  Figure  5  shows  the  measured  size  distribution  of  the 
sample  with  a  diameter  distribution  at  around  0.06  pm  to  0.125  pm.  Figure  6  shows  the 
DSC  curve  of  the  sample.  A  slightly  broader  peak  at  around  20  °C  is  due  to  the  phase 
change  of  PCM.  It  can  be  observed  that  the  phase  change  peak  is  around  20  °C  and  not 
between  28  °C  to  30  C,  typical  to  octadecane  which  is  probably  induced  by  the  small 
size,  where  the  surface  effect  is  much  more  significant  than  that  at  micro  scale. 


The  proposed  working  fluid  is  a  kind  of  two  component  slurry  that  consists  of  the 
MEPCMs  mixed  homogeneously  in  a  conventional  sensible  heat  transfer  fluid.  The 
assumption  of  homogeneous  Newtonian  fluid  is  valid  for  the  slurry  flow  mechanisms  for 
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the  particle  loading  range  of  0-25%  volume  [8].  The  effective  specific  heat  during  phase 
change  given  by  [6]: 


(1) 


hulk 


c 

where  p  is  the  weighted  average  specific  heat  given  by 


fraction  of  particles  undergoing  phase  change,  c  is  the  loading  fraction,  AT  is  the 
temperature  rise  of  fluid,  hSf  is  the  latent  heat  of  melting  of  PCM,  cpwf  is  the  specific  heat 
of  working  fluid  and  Cpcm  is  the  specific  heat  of  solid  PCM.  As  an  example,  assuming  the 
PCM  has  completely  change  its  phase  from  solid  to  liquid  at  the  exit  of  channels  (x=l), 
working  fluid  as  water,  the  loading  fraction  is  0.2  (c  =  0.2),  latent  heat  of  PCM  as 
250kJ/kg,  Cpcm  as  2  kJ/kgK  and  AT  as  3°C,  the  effective  specific  heat  will  be  21  kJ/kg.K. 
This  implies  that  the  specific  heat  is  increased  by  about  5  times  with  the  presence  of 


MEPCM. 


Figure  4:  Light  scattering  of  nanoPCM  in  right  jar;  left  jar  has  water 
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Figure  5:  Size  distribution  of  NEPCM particles  in  slurry  sample 
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Figure  6:  DSC  curve  of  NEPCM  slurry 

Colvin  and  Mulligan[9]  have  reported  that  with  increase  of  specific  heat  of  fluid  by  PCM 
slurry,  for  the  same  temperature  rises,  the  flow  rate  can  be  reduced  to  one-tenth  (or  by 
90%)  and  the  pumping  power  to  one-hundredth  (or  by  99%)  for  a  given  heat  flux.  Similar 
reductions  in  pumping  power  have  also  been  reported  by  Goel  et  al.  [10].  MEPCM  can 
greatly  enhance  the  performance  of  a  microchannel  heat  sink  by  providing  higher 
convective  heat  transfer  coefficients,  better  axial  temperature  uniformity,  and  reduced 
coolant  flow  rate  requirements. 


We  proposed  to  use  MEPCM  slurry  as  the  heat  transfer  fluid  inside  the  microchannels 
instead  of  only  liquid  for  the  following  advantages: 

.  High  rate  of  heat  transfer  per  unit  volume  due  to  the  ratio  of  surface  area  to 
volume  of  microchannels  leads  to  capability  of  removal  of  high  heat  flux. 

•  The  presence  of  MEPCM  particles  in  the  fluid  increases  the  overall  heat 
capacity  and  hence  the  mass  flow  rate  can  be  significantly  decreased  while 
maintaining  a  near  isothermal  condition  at  the  heat  source. 

•  The  reduction  in  mass  flow  rate  would  decrease  the  overall  pressure  drop, 
resulting  in  a  small  pump  and  less  power  consumption  for  fluid  circulation. 

•  The  Microencapsulated  PCM  shows  a  higher  thermal  cycling  resistance. 
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3.1a  Related  Work 


The  use  of  phase  change  materials  in  liquids  has  been  studied  extensively  for 
performance  enhancement  in  conventional  size  channels  [8,  12].  Colvin  et  al.  [11] 
reported  specific  heat  increase  up  to  five  times  and  heat  transfer  coefficient  increase  of  up 
to  2.8  times  in  flows  of  MEPCM  suspension.  Increase  in  convective  heat  transfer 
coefficient  partially  comes  from  the  increase  in  effective  specific  heat  due  to  the  latent 
heat  of  the  MEPCM.  Since,  Nu  oc  Pr04,  the  convective  heat  transfer  coefficient  for  a 
single-phase  fluid  can  be  described  in  terms  of  thermal  conductivity  and  specific  heat  by 
following  relation: 

h  oc  k06c°p4  (2) 

If  the  latent  heat  is  viewed  as  a  form  of  specific  heat,  the  increase  in  effective  specific 
heat,  which  includes  the  effect  of  latent  heat,  can  lead  to  an  increase  in  the  heat  transfer 
coefficient.  The  effective  thermal  conductivity  of  static  dilute  suspension  can  be 
evaluated  from  Maxwell's  relation  [12]: 

2  +  k„l  kf  +  2 c(k.  I  k  f  - 1) 

kh  = - p- - - -  (3) 

2  +  kp  /  kf  -c{kp  /  kf  - 1) 

Sohn  and  Chen  [13]  showed  that  the  thermal  conductivity  of  solid-liquid  suspensions 
increased  effectively  in  a  laminar  flow'  because  of  the  effects  of  micro-convection  around 
solid  particles  and  particle-to-particle  interaction.  It  was  found  that  the  degree  of  the 
enhancement  of  thermal  conductivity  increases  as  particle  diameter  increases.  For 
particles  of  diameter  smaller  than  100  pm,  this  enhancement  can  be  negligible.  An 
associated  heat  transfer  enhancement  is  therefore  believed  to  be  marginal  compared  to 
that  related  to  latent  heat  [14], 


The  primary  parameter  that  affects  the  heat  transfer  enhancement  of  the  phase  change 
slurry  flow  is  the  bulk  Stefan  number  [8].  It  is  defined  as  the  ratio  of  the  sensible  heat 
capacity  of  the  suspension  to  its  latent  heat  capacity.  For  the  constant  heat  flux  condition, 
it  can  be  expressed  as  [8]: 
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(4) 


Steh 


CpeiaSjlh) 

cL(Pplph) 


where  qw  is  the  wall  heat  flux,  ra  is  the  channel  radius,  pp  and  pb  are  the  density  of  the 
particle  and  the  bulk  fluid.  Experimental  investigation  demonstrated  [10]  that,  the  wall 
temperature  rise  decreases  with  the  decrease  in  Stefan  number  for  the  boundary  condition 
of  constant  wall  heat  flux.  From  the  definition  of  the  Stefan  number,  it  is  clear  that  for 
given  fluid  and  phase  change  materials,  its  value  depends  on  the  channel  size  and 
concentration.  Thus  for  best  cooling  results,  a  low  Stefan  number  should  be  maintained 
as  far  as  possible  by  either  by  increasing  the  microcapsule  concentration  or  by  reducing 
the  duct  radius. 


There  are  very  limited  experimental  investigations  on  heat  transfer  of  MEPCM  slurry. 
Goel  et  al.  [10]  conducted  an  experimental  study  for  cases  of  flow  in  a  circular  tube  of 
3.14mm  inner  diameter  with  a  constant  heat  flux  boundary  condition.  Yamagishi  et 
al.f  1 5]  experimentally  investigated  the  hydrodynamics  and  heat  transfer  characteristics  of 
slurry  containing  microencapsulated  phase  change  material  for  use  as  a  heat  transfer  fluid 
in  a  circular  tube  of  10mm  with  uniform  heat  flux  condition.  Goel  et  al's  experiment  is 
the  smallest  tube  size  we  can  found  up  to  now. 


As  mentioned  above,  small  tube  size  can  improve  cooling  performance  for  PCM  slurry 
flow.  Tao  et  al.  [16]  proposed  a  heat  sink  design  that  contains  a  3D  micro/nano  network 
and  liquid  fluid  mixed  with  nanosize  phase  change  materials  (NPCMs).  Hao  and  Tao  [17] 
developed  a  numerical  model  to  investigate  the  heat  transfer  characteristics  of  PCM 
slurry'  in  microchannels.  The  model  was  based  on  the  continuum  model  for  both  the 
carrier  fluid  and  PCM  particles.  The  interaction  of  momentum  between  liquid  and  solid  is 
coupled  based  on  the  empirical  correlation  for  packed  beds  or  fluidized  beds,  depending 
on  the  concentration  of  the  slurry.  The  authors  assumed  local  thermal  equilibrium  for 
both  liquid  and  solid.  Separate  volume-averaged  energy  equations  were  written  for  the 
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liquid  and  PCM  particles.  Two  energy  equations  are  coupled  by  an  interfacial  heat 
transfer  coefficient  estimated  based  on  Wakao  and  Kaguei's  relation.  They  compared 
their  simulation  [18]  with  experimental  result  for  a  microencapsulated  PCM  suspension 
flow  in  a  3mm  diameter  tube  [10]. 


In  order  to  reduce  the  temperature  rise  for  removal  of  high  heat  fluxes,  high  loading  is 
expected  to  increase  the  overall  specific  heat  of  in  working  fluid.  Hydrodynamic  and  heat 
transfer  characteristics  of  microencapsulated  PCM  slurry  in  microchannels  will  be  quite 
different  from  those  of  dilute  slurry  in  conventional  size  channels.  There  are  some 
fundamental  issues  to  be  addressed  for  understanding  flow  and  heat  transfer  of  PCM 
slurry  in  the  microchannel  as  a  first  step  towards  application.  The  numerical  and 
experimental  work  that  has  been  performed  till  now  is  described  in  the  next  two  sections. 
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3.2  Particle  Migration  in  Suspension  Flows 


Last  year,  the  formulated  LBM  code  was  used  to  simulate  the  particle  paths  when  the 
duct  shape  has  different  aspect  ratios  and  few1  results  were  presented.  This  year,  the 
simulations  were  continued  and  more  results  are  presented  after  a  brief  introduction  to 
particle  migration  and  LBM  methodology. 


Most  of  the  numerical  models  simulating  MEPCM  slurry  flows  assumed  homogeneous 
distribution  of  particles.  The  model  proposed  by  Tao  et  al.  [16]  use  fluidized  bed 
correlations  to  solve  the  particle  distribution  and  is  derived  only  for  two-dimensional 
flows.  So  it  is  important  to  understand  the  behavior  of  particles  in  microchannels  where 
the  length  of  the  channel  is  not  enough  for  the  particle  distribution  to  be  fully  developed. 


When  a  slurry  or  suspension  flow  inside  a  channel,  it  was  observed  by  many  researchers 
that  the  particles  tend  to  migrate  and  the  flow  becomes  nonhomogeneous  even  though  it 
is  homogeneous  at  the  inlet.  This  effect  has  been  first  observed  by  Siegre  and  Silberberg 
[19]  and  noted  that  a  neutrally  buoyant  particle  tends  to  migrate  to  an  equilibrium 
position  approximately  at  0.6r,  at  low  Re,  where  r  is  the  pipe  radius.  This  effect  is  known 
as  Siegre  and  Silberberg  effect  or  tubular  pinch  effect  (figure  7).  The  mechanisms 
suggested  for  particle  migration  in  concentrated  suspensions  include  hydrodynamic 
interactions,  electrostatic  interaction,  and  other  surface  interactions  that  become 
important  as  the  particles  are  close. 
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Flow 


Tube  Inlet  During  Migration  Ful|y  Developed  Profile 

Figure  7:  Particle  distribution  over  a  cross  section  of  tube 

There  has  been  growing  interest  in  the  flow  of  solid  suspension  and  the  particle 
distribution  recently.  In  cases  where  Brownian  motion  is  negligible,  the  migration  of 
particles  has  been  attributed  to  the  presence  of  inertial  forces.  However,  under 
vanishingly  small  particle  Reynolds  numbers,  inertial  forces  are  not  dominant  and  the 
migration  of  particles  depends  purely  on  the  shear-induced  forces.  Most  studies  on 
particle  migration  in  semi-concentrated  or  concentrated  suspension  have  been  focused  on 
the  cases  of  vanishingly  small  particle  Reynolds  number,  Rep  [20-30].  These  models  are 
based  on  different  simulation  techniques  such  as  Stokesian  dynamics,  Lattice  Boltzmann 
methods,  dissipative  particle  dynamics  and  the  Lagrange  multiplier  fictitious  domain 
method. 

3.2a  Simulation  of  Particle  Distribution  Using  LBM 

Lattice  Boltzmann  (LB)  methods  belong  to  the  family  of  mesoscopic  methods.  Each 
conservation  law  is  related  to  a  microscopic  quantity  which  is  conserved  exactly  by  the 
collision  operator  of  the  evolution  equation.  The  evolution  equation  describes  the 
dynamics  of  distribution  functions  moving  with  discretized  velocities  between  the  nodes 
of  the  computational  grid.  The  method  proposed  by  Ladd  [31,  32]  for  simulation  of 
particulate  suspensions  using  lattice  Boltzmann  method  (LBM)  has  been  used  in  the 
current  project  to  simulate  the  particle  distribution.  One  of  the  important  features  of  the 
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present  method  is  that  the  computational  cost  scales  linearly  with  the  number  of  particles 
in  the  domain.  In  the  lattice  Boltzmann  method  of  simulating  such  particulate 
suspensions,  the  discrete-velocity  models  of  the  fluid  combine  the  features  of  a  fully 
molecular  simulation  of  both  solid  and  liquid  phases,  at  the  same  time  being  orders  of 
magnitude  faster  than  computationally  intensive  solutions  of  the  Navier-Stokes  equations. 
As  compared  to  Navier-Stokes  based  solvers,  lattice  Boltzmann  simulations  are  fast, 
flexible  and  simple. 


The  Boltzmann’s  kinetic  equation  is  a  well  established  mathematical  model  of  a  fluid  at 
the  microscopic  level  which  describes  the  evolution  of  the  single  particle  distribution 
function.  The  problem  is  to  determine  the  time  evolution  of  the  particles  with  velocity  v 
around  the  space-time  point(r,/)  given  as/(r,v,/).  The  one-particle  distribution 

function  /  (r,/) describes  the  number  density  of  particles  at  a  particular  node  of  the 
lattice  rand  time  t  with  velocities  given  bye,.  The  equation  for  the  particle  distribution 
function  using  the  lattice  BGK.  (LBGK)  equation  is  given  by: 


f){\+el8t,t+ St)- 


(5) 


where  f‘q  is  the  equilibrium  distribution  of  particles  moving  in  direction  7’. 

The  density  per  node  and  the  macroscopic  momentum  flux  are  defined  in  terms  of  the 
particle  distribution  functions  by: 

(6) 


The  details  of  the  model  can  be  found  in  [33,  34],  LBM  was  used  to  simulate  distribution 
of  particles  in  a  channel  over  a  period  of  time. 
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3.2b.  Results  Using  LBM 


In  the  current  simulations,  particle  size  is  varied  between  d=9.65x  and  d=16.45x  and  the 
duct  cross-section  varies  from  lOdxlOd  to  10dx40d  for  a  range  of  aspect  ratios 
betw  een  1  and  4.  The  length  of  the  duct  is  kept  constant  at  20 d  for  all  the  simulations.  As 
the  lattice  Boltzmann  equation  is  weakly  compressible,  simulations  at  higher  Reynolds 
number  were  conducted  with  particles  of  size  d=16.45x  and  with  low  viscosity  fluids. 
Secondary  flow  is  know  n  to  emanate  from  the  comers  of  the  duct  and  can  propagate  into 
the  main  flow  [35].  The  magnitude  of  velocities  in  the  region  of  secondary  flows  is 
proportional  to  the  velocity  of  the  flow.  At  high  Re,  to  minimize  the  effect  of  the 
secondary  flows  simulations  were  conducted  with  low  viscosities.  The  Reynolds  number 
used  in  the  current  simulations  is  based  on  the  shorter  side  of  the  channel  cross-section, 
and  is  defined  as  Re  =  pHUavgj fu  . 

Results  in  square  microchannels  (aspect  ratio  =  1) 

Simulations  were  conducted  to  validate  the  results  for  a  square  duct  with  16-20  particles 
distributed  randomly  in  the  channel.  Figure  8  compares  the  randomly  assigned  initial 
particle  positions  with  the  positions  after  migration  was  found  to  be  completed,  around 
t=10,000ts  for  Re=100.  Here  ts  is  the  Stokes  time  and  is  given  by  ts=d/2UaVg.  It  can  be 
observed  that  the  particles  migrate  to  a  region  shown  by  a  dotted  circle  of  radius  given  by 
R=0.38,  where  the  non-dimensionalization  of  distances  is  done  based  on  the  shortest  side 
of  the  channel  cross-section,  from  the  center  of  the  square  duct  (figure  8(a)).  The  final 
particle  positions  are  found  to  be  at  an  average  distance  of  0.12  away  from  the  wall, 
which  is  consistent  w  ith  the  observations  of  Chun  and  Ladd  [35],  who  report  the  corners 
and  region  near  the  centers  at  a  distance  of  0.13  from  the  walls  to  be  stable  equilibrium 
positions.  In  figure  8(b),  the  particle  trajectories  for  Re=100  are  shown.  It  can  be 
observed  that  the  particles  near  the  walls  are  ‘repelled’  away  from  it,  due  to  a  lubrication 
force  that  is  generated  because  of  the  viscous  forces  generated  in  the  thin  film 
sandwiched  between  the  particles  and  the  wall.  On  the  other  hand,  particles  close  to  the 
center  have  a  constant  outward  movement  directed  away  from  the  center  of  the  channel. 
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This  inward  movement  can  be  explained  based  on  the  curvature  of  the  velocity  profile, 
which  changes  rather  smoothly  near  the  center  and  is  sharp  away  from  the  center. 


Y  Y 

(a)  (b) 


Figure  8:  Final  equilibrium  positions  at  10,000ts  for  Re=100  and  aspect  ratio  1 

channel 

Simulation  of  neutrally  buoyant  particles  at  Re=250  for  the  square  cross-section  channel 
gives  a  slightly  different  picture.  Figure  9  shows  such  a  case.  Particles  were  observed  to 
move  to  the  comers  or  a  region  close  to  the  center  of  the  sides.  As  compared  to  Re=100 
(figure  8),  the  particles  were  observed  to  be  more  localized  at  certain  preferred  zones  in 
the  channel.  More  so,  particles  were  found  to  be  aligned  in  chains  along  the  direction  of 
the  flow.  Figure  10  shows  the  three-dimensional  view  of  the  channel  and  the  particle 
locations.  Chains  of  particles  were  found  to  be  formed  at  an  average  distance  of  0.107 
from  the  walls  of  the  channel.  This  observation  is  qualitatively  consistent  with  the 
simulation  results  of  Chun  and  Ladd  [35],  who  observed  a  similar  behavior  at  Re=300. 
Their  results  show  that  the  length  of  the  stable  region  along  the  sides  of  the  channel 
decreases  as  Re  is  increased  from  300  to  500.  For  Re=500,  particles  were  shown  to  form 
uneven  spaced  clusters,  and  reside  near  the  corners  of  the  square  cross-sectioned  channel. 
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Figure  9:  Equilibrium  locations  of  particles  at  Re=250  and  aspect  ratio  l  channel 
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Figure  10:  Formation  of  chains  along  the  direction  of  the  flow  for  aspect  ratio  l  and 

Re-250 
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Results  in  microchannels  of  aspect  ratio=2 


Simulation  results  of  inertial  migration  of  neutrally  buoyant  spheres  at  an  aspect  ratio  of  2 
are  described  in  this  section.  Sixteen  particles  of  size  d/H=  1/10  were  distributed 
randomly  in  the  channel  initially.  Figure  1 1  represents  the  initial  and  final  positions  of 
these  particles  for  Re=100  and  compares  these  positions  once  equilibrium  is  reached.  It 
can  be  observed  that  the  particles  tend  to  be  spread  over  the  shorter  and  longer  sides  of 
the  channel,  and  are  not  localized  to  a  certain  point  or  region.  The  particles  align  along 
the  longer  side  at  an  average  distance  of  (0.176)  from  the  wall.  Along  the  shorter  side,  the 
particles  were  observed  to  be  positioned  at  an  average  distance  of  (0.171)  from  the  wall. 
Compared  with  the  equilibrium  locations  obtained  for  a  square  cross-sectioned  channel 
for  Re=l00.  it  is  apparent  that  the  particles  are  pushed  farther  away  from  the  wall  as  the 
aspect  ratio  is  increased. 
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Figure  11:  Particle  locations  for  Re=100  and  aspect  ratio  2 


Simulation  results  for  a  slightly  higher  Reynolds  number  of  250  show  a  slightly  different 
trend.  The  migration  is  faster  as  compared  to  Re=100.  as  the  particles  experience  a  higher 
fluid  inertia,  which  causes  them  to  move  faster  in  the  lateral  direction.  Also,  the  number 


19 


density  of  particles  near  the  comers  of  the  channels  is  found  to  increase  as  compared  to 
the  former  case.  Moreover,  particles  tend  to  move  away  from  the  longer  side  of  the 
channel.  This  indicates  the  presence  of  an  inner  ring,  as  shown  in  figure  12,  which 
becomes  apparent  for  the  first  time.  In  a  recent  experimental  study  [36],  the  presence  of 
an  inner  ring  of  particles  for  a  circular  channel  was  reported  for  Re>600  for  D/d=9-17.  In 
figure  12,  it  can  be  observed  that  the  particles  migrate  to  two  equilibrium  positions,  given 
by  0.14  and  0.25  from  the  walls  at  5600ts.  The  outer  ring  moves  more  close  to  the  wall 
for  Re=250  as  compared  to  Re=100,  consistent  with  the  arguments  laid  by  the 
perturbation  theories. 
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Figure  12:  Particle  locations  for  Re=250  and  aspect  ratio  2 

As  the  Reynolds  number  is  increased  further,  it  was  observed  that  the  region  around  the 
centers  of  the  longer  faces  is  depleted  of  any  particles  (see  figure  13).  Instead  these 
particles  move  inwards  towards  the  center  of  the  channel,  thereby  pushing  the  inner  ring 
farther  away  from  the  walls.  The  outer  ring  of  particles  is  found  to  be  at  a  distance  of  0.27 
from  the  wall.  On  the  other  hand,  particles  near  the  vertical  walls  tend  to  move  outwards 
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and  ultimately  reside  at  a  location  of  0.1  from  the  walls.  The  particle  concentration  was 
also  found  to  be  increase  near  the  comers. 


O  t  =  0 


Figure  13:  Particle  locations  for  Re=500  and  aspect  ratio  of  2 


Results  in  microchannels  of  aspect  ratio  =  4 

Simulation  results  of  inertial  migration  of  neutrally  buoyant  spheres  in  a  channel  of 
aspect  ratio  4  are  described  in  this  section.  Similar  to  Re=250,  16  particles  of  size 
d/H=  1/10  were  distributed  randomly  in  the  channel  initially.  Different  pressure  gradients 
were  employed  to  drive  the  Poiseuille  flow.  Figure  14  represents  the  migration  pattern  of 
these  neutrally  buoyant  particles  for  Re=100.  As  can  be  observed,  particles  near  the  walls 
are  pushed  inwards  and  particles  near  the  center  move  outwards.  At  equilibrium,  particles 
were  found  to  reside  at  a  distance  of  0.2  from  the  walls.  This  distance  is  more  as 
compared  to  the  same  inertia  case  for  aspect  ratio  of  2,  where  the  particles  were  found  to 
be  at  a  distance  of  0.17,  indicating  the  role  of  wall  forces  exerted  from  the  longer 
dimension. 
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At  a  slightly  higher  Reynolds  number  of  250,  a  behavior  similar  to  the  same  Re  case  of 
A=2  is  observed  (see  figure  15).  Particles  were  found  to  be  repelled  away  from  the  longer 
side  wall.  A  fraction  of  the  particles  were  observed  to  migrate  to  positions  that  were  0.13 
away  from  the  walls,  forming  an  outer  ring.  The  rest  of  the  particles  reside  at  a  distance 
of  0.27  from  the  walls.  A  quantitative  comparison  of  the  same  case  for  A=2  reveals  that 
both  the  inner  and  outer  rings  are  not  at  the  same  locations  for  the  two  different  aspect 
ratio  cases.  For  A=2,  the  outer  (inner)  ring  was  found  to  be  at  a  distance  0.14  (0.25)  from 
the  walls,  whereas  for  A=4,  this  distance  was  calculated  to  be  0.13  (0.27). 


To  compare  this  shifting  of  the  inner  and  outer  rings  with  change  in  Re  and  aspect  ratio, 
results  of  a  simulation  conducted  for  Re=500  is  shown  in  figure  16.  For  this  flow,  the 
number  density  of  particles  in  the  inner  ring  is  higher  than  the  number  density  of  particles 
that  migrate  to  the  walls  of  the  channel.  The  presence  of  an  inner  and  outer  ring  is  also 
clear,  with  these  regions  located  at  a  distance  of  0.3  and  0.1  from  the  walls  respectively. 


O  t  =  3363  ts 
O  t  =  6500  ts 


Figure  14:  Particle  positions  at  different  time  instants  for  Re-100  and  aspect  ratio  of  4. 
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Figure  15:  Particle  positions  for  Re=250  and  aspect  ratio  of  4 
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Figure  16:  Particle  positions  at  different  times  for  Re=500  and  aspect  ratio  of  4  -  The 
horizontal  dotted  lines  indicate  the  equilibrium  locations  along  which  the  particles 
seem  to  align,  at  0. 1  and  0.9  (outer  ring)  and  0.3  and  0. 7  (inner  ring) 


Experiments  have  revealed  that  as  Re  is  increased,  the  probability  of  finding  a  particle  in 
the  inner  ring  is  higher  as  compared  to  the  outer  ring  [36].  More  so,  for  D/d=42,  the 
presence  of  the  inner  ring  was  not  observed  until  Re>1200,  and  the  inner  ring  was  found 
to  be  centered  at  a  distance  of  about  half  of  the  radius  from  the  center  of  the  circular 
channel.  Current  simulation  results  indicate  a  similar  behavior  for  H/d=10.  For  Re=500 
and  A=4,  the  inner  ring  is  found  to  be  located  at  a  distance  of  0.3  from  the  walls  (0.2 
from  the  center).  In  addition,  the  number  of  particles  that  migrate  to  the  inner  ring  is  more 
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than  twice  the  number  of  particles  in  the  outer  ring.  With  an  increase  in  inertia  of  the 
flow,  the  inner  ring  is  found  to  move  towards  the  center  of  the  channel  between  Re=250 
and  Re=500. 


Asymptotic  theory5  predicts  that  for  planar  Poiseuille  flows,  the  particles  move  to  a 
location  that  is  0.12,  0.09  and  0.08  far  from  the  walls  for  Re=100,  250  and  500 
respectively.  Table  1  summarizes  the  results  obtained  using  LBM  simulations  for 
different  aspect  ratios  with  the  theoretical  analysis  for  planar  Poiseuille  flow  [37].  To 
summarize,  particle  migration  locations  are  dependent  on  the  aspect  ratio  and  the  inertia 
of  the  flow.  As  the  Reynolds  number  is  increased,  the  particles  move  to  an  equilibrium 
location  that  is  closer  to  the  wall.  On  the  other  hand,  as  the  aspect  ratio  is  increased,  the 
relative  equilibrium  distances  from  the  wall  increase  as  well.  For  Re>250  for  both  aspect 
ratios  of  2  and  4,  an  inner  ring  of  particles  was  formed.  The  location  of  this  inner  ring  is 
also  a  function  of  the  Reynolds  number  and  the  aspect  ratio.  With  an  increase  in  either  of 
the  two,  the  ring  was  observed  to  shift  to  a  location  that  was  farther  away  from  the  wall. 


Table  1:  Particle  equilibrium  locations  for  different  Re  and  aspect  ratios.  The  values 
of  the  inner  rings  are  mentioned  in  brackets. 


Re 

A=1 

A=2 

A=4 

Theoretical  Solution  [37] 

100 

0.12 

0.17 

0.2 

0.12 

250 

0.107 

0.14(0.25) 

0.13  (0.27) 

0.09 

500 

0.1  (0.27) 

0.1  (0.3) 

0.08 

3.2c  Conclusion 

In  this  work,  inertial  migration  of  neutrally  buoyant  particles  in  three-dimensional 
rectangular  channels  has  been  simulated  using  LBM.  Simulation  results  were  presented 
for  cross-section  aspect  ratios  of  1, 2  and  4,  for  100<Re<500  with  the  channel  to  particle 
size  ratio  fixed  at  10,  and  compared  with  existing  experimental  data  [19,  36]  asymptotic 
theory  [38],  and  simulation  results  for  a  square  channel  [35].  Migration  was  found  to  be 
faster  in  the  first  half  of  the  total  migration  time.  During  this  time,  the  particles  seem  to 
cut  the  flow  streamlines,  whereas  in  the  second  half  these  particles  move  parallel  to  the 


24 


walls  until  they  reach  a  final  equilibrium  position.  For  aspect  ratios  higher  than  1.  a 
difference  in  the  migration  location  was  observed.  For  a  fixed  Reynolds  number,  the 
distance  of  the  equilibrium  positions  from  the  walls  was  found  to  be  increase  as  the 
channel  aspect  ratio  was  increased.  For  Re=100.  all  the  particles  were  found  to  migrate  to 
an  outer  ring.  The  separation  between  this  particle  rich  zone  and  the  wall  was  found  to  be 
larger  as  compared  to  those  for  A=1 .  For  Rey  nolds  number  above  250,  the  presence  of  an 
inner  ring  was  evident,  whose  distance  from  the  wall  is  found  to  increase  with  increasing 
aspect  ratio.  For  a  fixed  aspect  ratio  and  Re>250,  the  particles  in  the  inner  ring  were 
found  to  be  pushed  further  away  from  the  wall  as  the  flow  Reynolds  number  was 
increased.  The  presence  of  the  inner  ring  has  previously  been  discussed  for  circular  and 
square  cross-sectioned  channels,  and  was  reported  to  appear  for  Re>500  with  H/d=10. 
Current  simulation  results  indicate  that  the  presence  of  an  inner  ring  is  possible  for 
channels  with  A>1  and  Re>250. 
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3.3.  Particle  Simulation  using  DFM 

As  mentioned  in  previous  section,  the  maximum  concentration  of  particles  that  was  used 
in  LBM  was  less  than  0.5%.  For  simulations  using  LBM  the  maximum  number  of 
elements  used  for  simulation  is  just  a  few.  Even  though  lattice  Boltzmann  methods  can 
provide  valuable  insight  to  many  body  interactions,  solutions  could  theoretically  be  found 
for  each  particle,  but  due  to  the  multibody  interactions  the  mathematics  becomes 
complicated  with  even  just  a  few.  The  number  of  particles  explicitly  modeled  in  these 
simulations  is  currently  limited  to  the  order  of  several  thousand  using  high-end  computers 
because  of  the  large  CPU  and  memory  requirements. 


The  particle  migration  also  depends  on  particle  concentration  and  presence  of  more 
number  of  particles  will  help  in  increasing  the  heat  capacity  of  the  fluid  during  melting. 
For  higher  concentrations,  DFM  was  used  to  solve  the  particle  distributions  and  the 
results  were  used  for  analyzing  the  thermal  performance  of  slurry.  This  model  was 
introduced  and  described  in  the  last  year  annual  report  and  the  results  obtained  this  year 
by  using  this  model  are  presented. 


3.3a  Diffusive  Flux  Model 

For  particle  suspensions  with  large  number  of  particles,  it  is  reasonable  to  consider  the 
suspension  as  a  continuum  for  most  applications  because  many  suspended  particles  are 
less  than  a  few  micrometers  in  diameter.  From  a  practical  point  of  view,  a  macroscopic 
constitutive  equation  is  preferable,  as  it  allows  the  modeling  of  realistic  macroscopic 
problems  that  contain  extremely  large  number  of  suspended  particles.  Several 
constitutive  models  have  been  put  forth  recently  which  fall  into  two  categories.  One 
category  comprises  of  models  based  on  conservation  of  mass  and  momentum  for 
suspension  components  [20],  the  other  category  comprises  models  based  on  shear- 
induced  particle  migration  and  diffusion  [26]. 
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The  second  category  models,  collectively  referred  as  diffusive  flux  model  (DFM)  are 
based  on  the  scaling  arguments  proposed  in  [24],  Leighton  and  Acrivos  [37]  suggested 
phenomenological  models  for  particle  migration  in  nonhomogeneous  shear  flow  typically 
due  to  spatial  variation  in  irreversible  interaction  frequency  and  effective  viscosity. 
Phillips  et  al.  [26]  adapted  the  scaling  arguments  of  Leighton  and  Acrivos,  and  proposed 
a  diffusive-flux  equation  to  describe  the  time  evolution  of  the  particle  concentration 
based  on  the  two-body  interactions. 


For  analyzing  the  particle  distribution  effect  on  thermal  performance  for  high 
concentrations,  the  Diffusive  Flux  model  (DFM)  proposed  by  Philips  et  al.  [26]  was  used. 

This  method  as  applied  to  a  three-dimensional  rectangular  duct  can  be  summarized  as 
follows.  The  shear-induced  migration  of  particles  is  the  result  of  shear  rate  gradients  and 
the  concentration  gradients.  The  particle  flux  based  on  the  phenomenological  model  is 
given  by: 


J  =  V2— V>7 

n 


(7) 


Kc  and  Kn  are  the  phenomenological  constants  that  must  be  determined  by  fitting  the 
predictions  of  the  model  and  the  experimental  results.  These  values  for  both  2-D  tube  and 
channels  flows  were  predicted  as  0.62  and  0.41.  The  dynamics  change  of  particle 
concentration  along  the  flow  is  governed  by  the  particle  balance: 

D<j> 


Dt 


=  V.J 


(8) 


With  the  assumptions  mentioned  above,  the  mass  and  momentum  equations  for  the  phase 
change  slurry  can  be  written  as: 

V.w  =  0  (9) 

V.[//(Vw  +  Vw' )]- Vp  =  0  (10) 

The  effective  viscosity  is  of  a  concentrated  suspension  can  be  given  by  the  Kriegefs 
formula  [27]: 
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n  =  v  f 


(ii) 


± 

tmj 


where  the  maximum  concentration  can  be  0.6  to  0.68  and  the  empirical  m  is  equal  to 
1 .82.  Diffusive  flux  model  is  valid  for  slow  flows  and  it  can  be  used  for  particle  Reynolds 
numbers  less  than  0.1  [28],  Kc  and  Kn  are  the  phenomenological  constants  that  must  be 
determined  by  fitting  the  predictions  of  the  model  and  the  experimental  results. 


Boundary  conditions  for  flow  and  particle  distribution 
For  particle  distribution,  the  follow  ing  boundary  conditions  are  used 
(j>  =  <j>Q  ,  at  the  inlet  (12) 

[tj(Vu  +  Vu'  )]n  =  0  p  =  p0  ,  particle  outlet  (13) 

1  3  71 

n.[K  c<f>V  [y<f>)  +  K  ( y<j> 2) - -V tf>]  =  0,  the  total  flux  at  the  wall  and  the  symmetric 

n  $0 

boundary  is  zero.  (14) 


Developing  length  of  particle  distribution 

Particle  distribution  can  be  considered  to  achieve  the  equilibrium  position  if  there  is  no 
change  in  the  particle  profile.  Nott  and  Brady  [20],  analyzed  the  time  scale  for  particles  to 
reach  steady  state  based  on  the  shear-induced-diffusion  hypothesis  [39]  and  is  given  by: 

t 

where  W  is  half  of  the  channel  width,  a  is  the  particle  radius,  <u>  is  the  average  velocity. 
For  dense  suspensions,  (<p>0.3),  12d(cp)~l.  Hence  the  length  along  the  channel  required 
to  achieve  steady  state  is  given  by: 


2  W 


a 


a  J  \2d(<f>)  <  u  > 


(15) 
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(16) 


_ 

2 W  \2d(</>)\  a  ) 


This  length  is  the  characteristic  length  scale  for  the  process  and  it  requires  several 
transition  lengths  before  the  final  equilibrium  state  is  reached.  Experiments  by  Hampton 
et  al.  [29]  shows  a  range  of  values  for  the  exponent  n,  increasing  linearly  from  0.4  for  cp  = 
0.2  to  1.8  for  cp  =  0.45.  Authors  in  [30]  found  that  the  value  of  n  is  around  1  for  particle 
concentration  of  0.24  and  around  2  for  particle  concentration  of  0.35.  Assuming  H  =  100 
pm,  a  =  2.5  pm  and  particle  concentration  is  0.3,  which  is  between  0.24  and  0.35,  particle 
steady  state  profile  can  be  achieved  between  a  length  of  4  cm  and  16  cm.  It  can  be 
observed  the  migration  development  length  is  very  long  compared  to  the  usual  length  of 
MMC  heat  sinks,  which  is  around  1mm.  Hence,  the  effect  of  particle  migration  on 
thermal  performance  was  investigated  in  traditional  microchannel  heat  sinks. 


For  the  present  study,  the  particle  distribution  is  modeled  in  microchannels  of  large 
aspect  ratio  (>=10),  where  the  flow  can  be  considered  to  be  two-dimensional.  This  was 
done  in  order  to  use  the  readily  available  2-D  values  for  Kt  and  Kn.  This  is  justifiable 
since  the  main  aim  of  this  work  is  to  look  at  the  thermal  performance  of 
nonhomogeneous  slurry  in  a  microchannel. 


Assumptions  for  thermal  simulation 

For  solving  the  temperature  of  slurry,  the  problem  was  formulated  based  on  the  follow  ing 
assumptions. 

i.  The  slurry  properties  are  a  function  of  particle  concentration.  Segre  and 
Silberberg  [19]  found  that  the  radial  migration  of  the  particles  is  a  function  of 
2.84th  power  of  particle/duct  diameter  ratio.  When  the  duct  or  channel  size  is  very 
small  like  in  the  case  of  microchannel.  the  ratio  can  be  large. 
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ii.  The  specific  heat  capacity  of  the  fluid  is  a  function  of  temperature.  The  slurry 
thermal  conductivity  is  a  function  of  the  particle  Peclet  number  and  varies  across 
the  flow  field. 

iii.  The  melting  inside  EPCM  particles  takes  place  over  a  range  of  temperatures, 
between  Ti  and  T2  with  the  peak  melting  point  at  a  temperature,  Tm. 

iv.  There  is  no  temperature  gradient  inside  the  particle  or  the  particle  melts 
instantaneously.  This  assumption  has  been  discussed  in  later  sections.  The  particle 
sizes,  where  this  assumption  is  valid  will  be  considered  for  simulation 

v.  The  effect  of  particle  depletion  layer  is  negligible.  The  particle  depletion  layer  is 
of  the  order  of  the  particle  radius  if  the  channel  size  to  particle  size  is  large  [40, 

41]- 

vi.  The  shape  of  the  encapsulated  particles  is  spherical.  The  shell  material  is  very  thin 
and  hence  the  particles  are  considered  to  consist  totally  of  the  phase  change 
material. 


Validation  of  assumptions  for  thermal  simulation 
•  Assumption  iv:  Particle  melts  instantaneously 

Since  the  length  of  the  channel  is  short  for  the  MMC  channel,  it  is  important  that  the 
PCM  particle  completely  melts  within  its  residence  time.  Charunyakorn  et  al.  [8]  applied 
the  method  proposed  by  Tao  [42]  to  calculate  solid-liquid  or  melt  interface  rp  (Figure  17) 
in  a  sphere  for  calculating  the  source  term  in  their  model. 


Solid  Wax 

Figure  1 7:  Particle  melting  process 
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Assuming  the  particle  has  to  melt  99%  by  the  time  it  exits  the  channel,  the  required 
temperature  difference  between  the  surrounding  fluid  Tf  and  the  melting  temperature  Tm, 
of  the  PCM  w  as  calculated  (Equation  17)  using  the  same  analogy. 


Tfm=Tf-Tm  = 
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(17) 


Where  the  Biot  number  of  the  particle  is  given  by: 


Bi„  =  - 
'  k 


iff 


2(1  -c) 


(18) 


pi'm  2-3 c3  +c 

Figures  18  and  19  show  the  required  temperature  difference  between  the  fluid  and  the 
particle  melting  temperature  with  water  and  PAO  as  base  fluids  in  a  channel  of  100  pm 
width  and  1  cm  length  for  different  Re.  It  can  be  observed  that  the  required  Tfm  increases 
with  the  increase  in  particle  radius.  It  can  be  observed  that  the  particle  may  not  be 
completely  melted  if  the  channel  length  is  short.  For  calculations  properties  in  Table  2 
were  used. 


Particle  Diameter  (micron) 


Figure  18:  Tf-  Tmfor  different  particle  diameters  (base  fluid  -  water) 
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Particle  Diameter  (micron) 

Figure  19:  T/—  Tmfor  different  particle  diameters  (base  fluid  -  PAO) 


Table  2:  Thermophysical  properties  used  in  Figures  3.2  and  3.3 


Density 

(kg/m3) 

Specific 

heat 

(J/kg.K.) 

Thermal 

conductivity 

(W/m.K.) 

Viscosity 

(kg/m.s) 

Latent 

heat 

(J/kg) 

PAO 

783 

2242 

0.143 

4.45x1  O'3 

- 

n-octadecane 

815 

2000 

0.18 

- 

244x1 03 

Water 

997 

4180 

0.604 

1  x  1  O'3 

- 

Governing  equations  for  thermal  simulations 


After  the  particle  distribution  is  solved  using  the  DFM,  the  thermal  performance  of  slurry 
will  be  modeled  assuming  slurry  as  a  bulk  fluid  with  varying  properties  which  are  a 
function  of  particle  concentration.  The  energy  equation  for  the  slurry  and  the  boundary 
conditions  are: 


f  dT  8T 
pc 3  u—  +  v— 
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8x  8y 


8x 


( 


kd-L 


\ 


8x )  8y 

T  =  Tjn  at  the  inlet  for  fluid 

q.n  =  (pCpUT).n;  convective  heat  flux  boundary  condition  at  the  outlet/exit 


(19) 

(20) 
(21) 
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q.n  =  qw;  constant  heat  flux  at  the  wall 
q.n  =  0;  insulation  at  the  symmetry  or  channel  half-width 


(22) 

(23) 


Effective  Thermal  Conductivity 

Slurry  thermal  conductivity  and  specific  heat  is  defined  follows.  For  dilute  suspensions  of 
static  bulk  fluids,  thermal  conductivity  of  the  suspension  can  be  defined  as  for 
conductivity  in  a  medium  with  distributed  spherical  particles  [12]: 


(24) 


Yamada  and  Takahashi  [43]  investigated  experimentally  the  thermal  conductivity  of 
suspensions  of  particles  of  different  shapes.  They  found  an  excellent  agreement  between 
their  experimental  values  for  suspensions  of  spherical  particles.  The  results  for  other 
shapes  compare  less  favorably.  For  flowing  slurries,  the  effective  thermal  conductivity  is 
higher  than  that  calculated  by  Equation  3.18  due  to  diffusion  related  enhancements.  For 
dilute  suspensions,  it  can  be  evaluated  as  follows  [8]: 


(25) 


/  =  I  +  BcPe"' 

1  p 

B  =  0,  m  =  1.5.  Pep<0.61 
fl  =  1.8  w  =  0.18  0.67  <  Pep  <  250 

B  =  3.0.  m  =  jy,  Pep  >  250 


The  particle  Peclet  number  is  defined  as, 


(26) 
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Since  the  velocity  is  not  fully  developed  in  the  current  analysis,  the  shear  rate  is  a 
function  of  all  the  spatial  coordinates  and  corresponding  velocities.  The  magnitude  of  the 
shear  rate,  e  can  be  calculated  using  the  following  equation: 


e  = 


^ZZ  r„r 

VZ  '  J 


J i 


J 


(27) 


y  is  the  shear  rate. 

Effective.  Density. 

For  slurries,  the  density  can  he  calculated  by  the  weighted  mean  method  given  below  [8]: 

Ph  =cPm  +0-c)P/  (28) 

The  mean  density,  pm  is  the  average  of  the  densities  of  PCM  solid  and  liquid  phases. 


Effective  Specific  Heat. 

The  phase  change  inside  the  particles  can  be  modeled  using  latent  heat  method  or 
effective  specific  heat  capacity  method.  Latent  heat  method  can  be  used  to  model  pure 
PCMs  where  as  the  effective  specific  heat  method  can  model  the  melting  range  of  PCMs 
easily.  Alisetti  and  Roy  [44]  have  used  various  profiles  for  the  specific  heat  of  PCM  for 
calculating  the  effective  specific  heat  and  have  shown  that  the  difference  between  the 
solutions  is  less  than  4%.  When  PCMs  are  encapsulated,  due  to  the  nucleation,  there 
usually  is  a  melting  range  for  the  PCMs,  and  this  increases  with  decrease  in  particle  size. 
Of  all  the  different  profiles,  sine  profile  (Figure  20)  has  been  used  to  represent  EPCM 
particle  specific  heat.  This  was  done  in  order  to  avoid  sudden  variation  in  the  property 
and  help  in  having  less  convergence  problems.  As  there  is  not  much  variation,  the 
specific  heat  of  solid  and  liquid  phases  of  PCM  is  assumed  to  be  equal. 
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The  effective  specific  heat  of  the  slurry  can  be  calculated  as: 


(29) 
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cp,b=cmcpp+{\-cm)cpf  (30) 

The  value  specific  heat  of  the  particle  in  Equation  30  is  equal  to  cppcm  when  the 
temperature  of  the  particle  is  outside  the  melting  range  and  is  given  by  Equation  29,  when 
the  particle  temperature  is  within  the  melting  range. 


Figure  20:  Specific  heat  of  EPCM,  function  of  temperature 


3.3b  Particle  Distribution  Results 

The  microchannel  dimensions  used  for  simulation  of  particle  distribution  are  100  pm 
width.  1mm  height  and  1  cm  length.  As  mentioned  before,  the  simulation  domain  can  be 
considered  2D  for  such  high  aspect  ratios.  Commercial  software  COMSOL  was  used  for 
all  the  numerical  simulations  in  the  work  [45].  Following  parameters  were  varied. 

i.  Base  fluid:  Two  different  fluids  were  used,  water  and  poly-alpha-olefin  (PAO). 
PAO  is  a  dielectric  fluid  used  for  cooling  of  military  avionics  applications.  It  is  a 
stable  and  inexpensive  fluid. 

ii.  Particle  diameter:  lOOnm,  1  pm,  5  pm 

iii.  Particle  concentration:  0.05,  0. 1 , 0.3 

iv.  Channel  Reynolds  number:  Since  DFM  is  valid  for  Rep<=0.1,  it  was  found  that 
the  maximum  Reynolds  number  of  the  fluid  that  can  be  modeled  was  calculated 
using  the  following  calculation: 
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Rep  =  Remax(a/H)2  (31) 

For  current  simulation,  a  =  2.5  pm,  H  =  100  pm.  Substituting  the  value  of  Rep  =  0.1, 
Remax  is  160.  Thus  three  different  channel  average  Reynolds  numbers  were  used  for 
simulating  the  particle  migration.  Table  3  below  shows  the  inlet  velocities  when  water 
and  PAO  were  used.  For  the  calculations,  the  density  and  viscosity  of  pure  fluid  are  used 
from  Table  2. 


Table  3:  Re  values  used  for  simulating  particle  migration 


Channel  Reavg 

^  in>  Water 

Vi„,  PAO 

80 

0.357  m/s 

2.276  m/s 

60 

0.268  m/s 

1 .707  m/s 

30 

0.134  m/s 

0.85  m/s 

Effect  particle  diameter  on  particle  migration 

Figure  21  shows  the  particle  distribution  at  a  length  of  1  cm  for  all  the  three  particle 
diameters  used  with  PAO  as  base  fluid.  The  particle  mass  concentration  is  0.3  and  Reavg 
is  80.  It  can  be  observed  that  the  migration  is  highest  for  5  pm  particle  and  there  is  no 
migration  for  100  nm  particle  size. 


Effect  of  particle  volume  concentration  on  particle  migration 

Figure  22  shows  the  particle  migration  at  a  length  of  1  cm  from  the  inlet  for  a  particle 
diameter  of  5  pm  and  Reavg  of  80.  It  can  be  observed  that  the  migration  is  high  when  the 
particle  mass  concentration  is  higher. 
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Nondimensional  distance  from  channel  center,  x/W 
Figure  21:  Particle  distribution  with  varying  particle  diameter  (PAO) 


Nondimensional  distance  from  channel  center,  x/W 

Figure  22:  Particle  distribution  as  a  function  of  concentration  (PAO) 


Effect  of  inlet  velocity  on  particle  migration 

Simulations  with  three  different  ReaVg  were  run  and  it  was  found  that  there  was  no  effect 
of  inlet  velocity  on  particle  migration.  Figure  23  shows  the  particle  concentration  profile 
obtained  with  water  based  slurry  for  all  three  Reavg.  This  type  of  result  was  also  obtained 
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in  [46].  The  reason  could  be  that  the  diffusive  flux  model  does  not  include  the  inertial 
forces  and  hence  the  variation  in  velocity  does  not  effect  the  particle  migration.  The 
variation  in  concentration  profile  was  only  dependent  on  the  a/H  ratio,  i.e.,  the  particle 
radius  to  the  tube  radius  or  channel  width. 


Nondimensional  distance  from  channel  center,  x/W 
Figure  23:  Particle  distribution  profile  (base  fluid  -  water) 


Particle  migration  along  the  length  of  the  channel 

It  is  interesting  to  observe  the  transformation  of  homogeneous  profile  into 
nonhomogeneous  along  the  length  of  the  channel.  Figure  24  shows  the  particle 
distributions  at  different  lengths  from  inlet  for  PAO.  The  Reavg  used  is  80  and  particle 
concentration  is  0.3.  It  can  be  observed  that  the  particle  migration  rate  increases  with 
increase  in  length. 
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NondimensionaJ  distance  from  channel  center,  x/W 


Figure  24:  Particle  concentration  at  different  locations  along  length  (PAO) 


3.3c  Thermal  Performance  Results 

The  effect  of  particle  migration  on  thermal  performance  was  investigated  by  comparing 
the  thermal  results  obtained  assuming  the  particle  profile  is  developing  (uniform  profile 
at  inlet  and  nonuniform  profile  at  the  exit  as  shown  in  Figure  25)  and  the  results  obtained 
assuming  the  particle  profile  is  homogenous  (uniform  at  both  inlet  and  exit  as  shown  in 
Figure  26). 


0 

0  0.2  0.4  0.6  0.8  1 


Nondimensional  distance  from  channel  center,  x/W 


Figure  25:  Particle  concentration  profile  for  developing  assumption 
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Figure  26:  Particle  concentration  profile  for  homogeneous  assumption 


The  particle  distribution  obtained  for  a  particle  diameter  of  5  pm  and  particle 
concentration  of  0.3  was  used  for  thermal  simulations  as  maximum  migration  was 
observed  at  these  conditions.  Figures  27  and  28  show  the  difference  in  maximum  wall 
temperature  and  the  inlet  temperature  for  both  fluids.  The  average  Re  used  was  80  and  a 
melting  range  of  6  K  was  used.  PCM  melting  peak  was  assumed  to  be  at  27  °C.  The 
difference  in  wall  temperature  for  both  assumptions  is  not  more  than  1.4  K  for  water  and 
is  1  K  for  PAO. 


Figure  2  7:  Difference  in  maximum  wall  temperature  and  inlet  temperature  for  both 

assumptions  (fluid  -  water) 
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90 


Figure  28:  Difference  in  maximum  wall  temperature  and  inlet  temperature  for  both 

assumptions  (fluid  -  PAO) 


Analysis  of  the  results 

From  previous  results,  it  can  be  understood  that  the  maximum  wall  temperature 
difference  for  both  particle  distribution  assumptions  was  not  more  than  1.4  K  for  all  the 
heat  fluxes  considered.  The  reason  can  be: 

•  Simultaneous  developing  of  both  concentration  and  temperature  profiles 


Since  the  concentration  profile  is  still  developing  along  the  length,  the  deviation  of 
concentration  from  the  average  value  (0.3)  is  not  immediate.  Figure  29  shows  the 
concentration  profile  at  a  distance  of  35  pm  from  the  center  (x/W  =  0.7)  along  the  length 
of  the  channel  from  inlet  for  PAO.  It  can  be  observed  that  the  minimum  value  of  particle 
concentration  at  the  wall  is  0.23,  which  happens  to  be  at  the  exit.  Figure  30  shows  the 
resultant  local  wall  temperature  from  inlet  to  exit,  assuming  both  homogeneous  and 
developing  particle  distributions.  It  can  be  observed  that  the  local  wall  temperature  for 
both  cases  is  same  up  to  a  certain  length  of  the  channel.  The  difference  in  temperature 
profiles  for  both  cases  will  be  the  combined  effect  of  thermal  boundary  layer  and  particle 
distribution  development,  which  is  negligible  for  a  certain  length  of  the  channel.  It  should 
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also  be  noted  that,  once  the  local  temperature  is  higher  than  the  melting  range,  PCM 
latent  heat  is  already  used  up  and  the  effect  of  migration  does  not  possibly  affect  the  wall 
temperature.  Hence,  the  difference  between  both  assumptions  did  not  affect  the  thermal 
performance  considerably. 


Figure  29:  Particle  concentration  profile  development  along  the  channel 
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Figure  30:  Local  wall  temperature  for  both  profile  assumptions  along  the  channel 
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Another  observation  from  the  results  is  that  the  maximum  wall  temperature  for 
developing  profile  assumption  is  lower  compared  to  maximum  wall  temperature  with 
homogeneous  assumption  in  case  of  water.  For  PAO,  this  behavior  was  found  to  be 
opposite.  The  reason  for  this  behavior  can  be  the  effective  thermal  conductivity  of  the 
fluid  at  the  wall.  In  case  of  water,  when  homogeneous  assumption  was  used,  the  thermal 
conductivity  was  lower  at  the  wall  resulting  in  less  heat  transfer  from  wall  to  fluid  and 
hence  high  wall  temperature  compared  to  developing  profile  assumption.  This  behavior  is 
reverse  for  PAO,  i.e.,  the  thermal  conductivity  is  higher  when  homogenous  assumption 
was  used  and  hence  lower  wall  temperature  compared  to  developing  profile  assumption. 


From  the  above  analysis,  it  can  be  observed  that  the  maximum  error  in  predicting  the  wall 
temperature  is  around  1 .4  K  for  the  conditions  used,  if  a  homogeneous  assumption  is  used 
instead  of  nonhomogeneous  assumption.  This  value  was  obtained  for  water  at  a  heat  flux 
of  20  W/cm2,  and  the  value  of  Twaii  -  Tm  is  29.2  K  for  developing  profile  and  30.6  K  for 
homogeneous  profile.  In  practical  applications,  this  small  temperature  difference  can  be 
neglected  owing  to  the  many  assumptions  used  in  numerical  simulations.  Unless  very 
accurate  prediction  is  needed,  including  the  particle  migration  effect  is  not  very  crucial. 


3.3d  Thermal  Performance  with  Fully  Developed  Profile  Assumption 

Since  a  major  difference  was  not  observed  in  developing  and  homogeneous  profiles, 
another  set  of  simulations  were  run  to  compare  the  thermal  performance  of  slurry  when 
the  particle  distribution  is  assumed  homogeneous  and  the  particle  distribution  is  fully 
developed  (Figure  31).  Particle  diameter  used  is  5  pm.  particle  concentration  is  0.3  and 
Remax  is  160.  Fluid  used  is  PAO.  For  parametric  study  the  effect  of  inlet  temperature, 
base  heat  flux,  melting  range  were  studied.  Figure  32  shows  the  concentration  profiles  at 
different  lengths  from  inlet.  From  the  figure,  it  can  be  observed  that  the  particle  migration 
rate  reduces  significantly  after  certain  length  from  the  inlet.  The  particle  concentration 
profile  varied  rapidly  from  inlet  up  to  2.5  cm  length  and  the  change  in  concentration 
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profile  variation  from  2.5  cm  to  10  cm  is  very  slow.  The  difference  between  profiles  at  a 
length  of  7.25  cm  from  inlet  and  10  cm  from  inlet  is  negligible.  Hence,  it  can  be  assumed 
that  the  particle  profile  is  fully  developed. 


Figure  31:  Particle  concentration  profile  for  fully  developed  assumption 
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Figure  32:  Particle  concentration  profiles  at  various  lengths  from  inlet  (base  fluid  - 

PAO) 
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Effect  of  fluid  inlet  temperature 


The  inlet  temperature  of  the  fluid  was  first  varied  at  a  heat  flux  of  20  W/cirf  and  melting 
range  of  5  K  (25.5  °C  to  30.5  °C,  and  peak  is  28  C).  Figure  33  shows  the  results.  It  can 
be  observed  that  an  inlet  temperature  of  26  C  (slightly  after  lower  melting  point)  will  be 
useful  in  obtaining  the  maximum  thermal  performance  at  any  profile  assumption.  The 
wall  temperature  difference  between  homogeneous  and  fully  developed  profile 
assumptions  is  around  3.6  K,  whereas  the  difference  between  homogeneous  and 
developing  profile  assumptions  is  1  K. 

Effect  of  wall  heat  flux 

Figure  34  shows  difference  between  maximum  wall  temperature  and  inlet  temperature  for 
all  the  three  assumptions.  It  can  be  observed  that  the  wall  temperature  difference  is 
maximum  between  fully  developed  and  homogeneous  assumptions  is  at  a  heat  flux  of  30 
W/crrT.  When  the  heat  flux  is  increased  to  40  W/cm2  and  50  W/cm\  the  bulk  temperature 
of  fluid  at  the  outlet  is  near  the  melting  end  temperature  (making  the  PCM  specific  heat 
compared  to  PAO)  and  hence  the  effect  of  migration  decreases. 
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Figure  33:  Tm,u.  mav  -  Tin  Tin  for  three  assumptions  at  different  inlet 

temperatures 
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Wall  Heat  Flux(W/cm2) 

Figure  34:  Twau  -  Tin  for  three  profile  assumptions  at  different  heat  fluxes 


Effect  of  PCM  melting  range 

Five  different  melting  ranges  were  considered  for  the  study.  Figure  35  shows  the  result  of 
varying  melting  range.  It  can  be  observed  that  thermal  performance  variation  between 
homogeneous  and  nonhomogeneous  assumptions  decrease  with  increase  in  melting 
range.  This  can  be  due  to  lowered  specific  heat  peak  (because  of  averaging  over  a  large 
melting  range)  which  makes  the  latent  heat  effect  of  PCM  comparable  to  specific  heat  of 
base  fluid.  Therefore,  the  effect  of  migration  of  PCM  particles  to  center  diminishes  with 
increase  in  melting  range. 


Since  a  typical  microchannel  can  also  be  2  cm  length,  few  more  simulations  were  run  to 
compare  the  thermal  performance  of  slurry  with  developing  profile  and  homogeneous 
profile.  Figure  36  shows  the  comparison  of  wall  temperature  for  both  assumptions  when 
the  channel  length  is  2  cm.  The  maximum  wall  temperature  difference  between  both 
profile  assumptions  was  1 .2  K  as  compared  to  1  K  when  the  channel  length  was  1  cm. 
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Figure  35:  Tw„ut  max  -  Ti„for  all  three  profiles  at  different  melting  ranges 

The  thermal  performance  of  slurry  was  also  analyzed  in  a  3D  channel  of  100  pm  width,  1 
mm  height  and  1  cm  length.  This  was  done  to  compare  the  maximum  wall  temperature 
obtained  in  3D  and  2D  simulations.  A  base  heat  flux  of  200  W/cm:  (equal  to  20  W/cm2, 
at  the  wall  for  2D)  was  used  for  3D  thermal  simulations.  Table  4  below  shows  the 
difference  between  2D  and  3D  results  for  all  the  three  assumptions.  It  was  found  that  the 
maximum  temperature  of  the  wall  (fluid-solid  boundary  in  width  direction)  is  slightly 
lower  than  the  wall  temperature  predicted  using  2D.  This  is  because  in  case  of  3D,  the 
increase  in  wall  heat  transfer  area  because  of  the  addition  of  the  base  results  in  a  lower 
wall  temperature. 
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Figure  36:  Results  at  a  channel  length  of  2  cm 


Table  4:  Thermal  results  obtained  in  3D  using  2D  particle  distribution 


Analysis 

T  wall,  11  NT 

dT  bulk,  IINT 

T wall,  Nil 

dTbulk,  NH 

Tw*u,  u 

dT  bulk,  n 

3D 

322.34 

300.9433 

324.561 

301.0379 

322.031 

300.8734 

2D 

323.409 

299.9423 

326.1888 

300.121 

322.491 

299.9458 

3.3e  Conclusion 

The  migration  was  found  to  depend  on  the  ratio  of  particle  diameter  to  channel  width. 
The  larger  the  particle  diameter,  the  more  rapid  is  the  particle  migration.  It  was  also 
found  that  the  particle  profile  varied  slightly  when  the  base  fluid  is  different.  It  was  found 
that  for  very  short  lengths  (~1  mm),  the  migration  of  particles  is  negligible. 


Thermal  simulations  using  the  results  of  diffusive  flux  model  showed  that,  when  the 
profile  is  developing,  the  maximum  wall  temperature.  Twaiimax  is  around  1.4  K.  different 
compared  to  Twaii,max  obtained  assuming  homogeneous  particle  distribution.  In  case  of 
water,  the  maximum  wall  temperature  assuming  developing  profile  is  lower  compared  to 
the  maximum  wall  temperature  obtained  assuming  homogeneous  distribution,  where  as 
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this  behav  ior  is  opposite  in  case  of  PAO.  The  reason  can  be  because  of  the  decrease  or 
increase  in  overall  thermal  conductivity  of  the  fluid  at  wall.  The  maximum  wall 
temperature  in  case  of  assuming  fully  developed  profile  showed  a  maximum  difference  of 
5  K  compared  to  that  of  assuming  homogeneous  distribution.  However,  for  obtaining 
such  profile,  the  length  of  channel  should  be  very  long  or  of  the  order  of  10  cm  in  the 
case  considered  which  is  not  likely  in  case  of  microchannels  (due  to  the  associated  high 
pressure  drops).  Hence,  it  can  be  concluded  that  it  is  not  important  to  include  the  particle 
migration  in  the  numerical  model. 
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3.4  Parametric  Study  in  MMC  Heat  Sinks 


The  parametric  study  results  presented  last  year  comprises  of  investigating  performance 
of  PAO  based  slurry  flow  in  MMC  heat  sinks.  This  section  discusses  the  effect  of  base 
fluid  thermal  conductivity  and  channel  dimensions  on  thermal  performance  of  slurry.  The 
model  is  briefly  introduced  here  and  parametric  study  results  are  presented  after 
comparing  with  results  obtained  in  heat  transfer  experiments.  From  section  3.3,  it  is 
evident  that  the  particle  migration  effect  is  negligible  and  hence  it  is  not  included  in  the 
model. 

3.4a  Numerical  Model 
Schematic  of  the  simulation  domain 

Figure  37  shows  the  schematic  of  the  heat  sink  domain  that  was  considered  for 
parametric  study.  The  flow  rate  at  the  entrance/inlet  of  the  channel  will  be  the  total  flow 
rate  at  the  heat  sink  inlet  divided  by  the  total  number  of  channels,  assuming  the  same 
amount  of  flow  rate  enters  each  channel.  Only  half  of  the  channel  was  considered  for 
simulation  due  to  symmetry. 


Figure  37:  Schematic  of flow  domain 
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Dimensions  used  are  shown  in  configuration  -  I  of  Table  5.  It  was  assumed  that  there  are 
441  microchannels  in  the  heat  sink  (the  specifications  of  MMC  heat  sink  used  for  heat 
transfer  experiments)  and  the  flow  rate  is  assumed  to  be  same  at  the  inlet  of  each  channel. 


Table  5:  Geometric  configurations  used  for  numerical  simulation  (units  in  jam) 


Configuration 

Base  fluid 

H 

2*tch 

L 

Hm 

2*tm 

Lm 

2X 

lbas« 

I 

Water,  PAO 

533 

101 

1000 

50 

201 

250 

352.52 

250 

11 

Water 

375 

25 

1000 

100 

100 

250 

75 

250 

Assumptions 

All  the  assumptions  mentioned  in  subsection  3.3  along  with  the  following  are  used  for 
formulating  the  problem. 

i.  The  particle  distribution  is  homogeneous.  This  is  valid  as  explained  in  Chapter  3. 

ii.  Particle  always  follow  the  fluid 

iii.  The  total  volumetric  flow  at  the  inlet  of  heat  sink  is  equally  divided  and  is  same  in 
each  channel. 

Validation  of  assumntions  major  assumptions 

•  Assumption  ii:  Particle  always  follow  the  fluid 

The  amount  of  time  taken  by  the  particle  to  reach  the  fluid  velocity  can  be  calculated 
using  Equation  32.  Using  the  drag  force  on  the  particle  (Figure  38) 

Fn=m^  =  6wfRr(ux-v)  (32) 


Figure  38:  Drag  on  a  spherical  particle  in  a  fluid  field 
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Integrating  the  equation  and  using  v  =  0  at  t  =  0,  to  find  the  integration  constant,  Equation 
32  becomes 


I  = 


2  PrRr~  Ini  K 


9  M, 


(««  ~v) 


(33) 


Defining  the  time  constant  as  in  Equation  4.3,  the  above  equation  becomes  Equation  35. 


r  =  • 


2  PrR 


r  p 


9  ft, 


It  -  V 


=  e 


(34) 

(35) 


At  t  =  t,  v  =  0.63  Uoo.  This  implies  that  the  particle  reaches  63  percent  of  fluid  velocity  in  i 
seconds.  The  time  constant  values  for  different  particle  diameters  with  water  and  PAO  as 
base  fluids  is  shown  in  Table  6.  For  the  values  below',  the  density  of  particle  used  is  774 
kg/m\  and  dynamic  viscosity  of  PAO  and  water  used  are  4.45x10°  Pa.s  and  7.39x1  O'4 
Pa.s  respectively. 


Table  6:  Time  constantfor  different  particle^  sizes 


Particle  diameter  (|im) 

^  water  (ns) 

x  pao  (ns) 

0.1 

0.57 

0.09 

1 

57.5 

9.64 

5 

1437.0 

241.0 

Governing  Equations 


With  the  assumptions  mentioned  above,  the  mass  and  momentum  equations  for  slurry  can 
be  written  as: 


du  dv  dw 
—  +  —  +  —  =  0 
dx  8y  dz 


(36) 


P 


dll  du  du  ' 

u —  +  v — +  vv — 

dx  dy  dz  j 


dp  d  (  du\ 
dx  dx  V  dx  J 


d  f  du ' 

+  —  // — 

dy{  dy) 


\ 

J 
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(37) 


(38) 


( 


dv  dv 


dv 


\ 


U - h  v - 1-  w — 

dx  dy  dz 


,_dp+d_ 


d 

M—\  + 


dy  dx{  dx  J  dy 


dv^ 

dy  j 


dz 


(  dw  dw  dw N 

P\  u - 1-  v - 1-  w — 

(  dx  dy  dz 


dp  d 
dz  dx 


cHv^l 

dx) 


+ — 

dy 


V  — 

dv 


+  ■ 


dz 


dv) 

V— 

.  dz) 


dw 


The  energy  equation  for  the  fluid  (slurry/pure  fluid)  is: 


( 


P°t 


dT  dT  dT 

u - K  v - 1-  rv — 

dx  dy  dz 


\ 


d_f 

~dx 


k?L 

dx 


+  - 


dy 


r  k— 
dy 


d  ( , dT^ 

+ - 

dz 


k  — 
dz 


Energy  equation  for  the  microchannel  wall/fin  is: 

d2T„,  d%  d%.  n 

■  + - T-  + - T-  =  0 


dx"  dy 


dz 1 


(39) 


(40) 


(41) 


The  corresponding  thermophysical  properties  of  pure  fluid  and  slurry  were  used  for 
Equations  36  to  40. 


Boundary  Conditions 

For  flow1  inside  the  channels: 
u  =  0  at  the  microchannel  and  manifold  walls, 

p  =  po,  atmospheric  pressure  at  the  outlet 

u  =  (0,  0.  |vt| )  at  the  inlet 

For  wall: 

q.n  =  qw;  constant  heat  flux  at  the  base  of  the  fin 
q.n  =  0;  adiabatic  on  all  other  outer  walls 
For  fluid: 

T  =  T,n  at  the  inlet 

q.n  =  (pcpuT).n;  convective  heat  flux  boundary  condition  at  the  outlet/exit. 


(42) 

(43) 

(44) 

(45) 

(46) 

(47) 

(48) 
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For  wall  and  liquid  interface: 


Tw  =  T;  continuity  of  temperature  (49) 

dT  dT, 

—  kw — —  =  —kf - ;  which  is  the  continuity  of  heat  flux  (50) 

dn  dn 

Dynamic  Viscosity ; 

For  suspensions  that  exhibit  Newtonian  characteristics,  the  bulk  viscosity  can  be 
calculated  from  an  empirical  formula  in  terms  of  the  particle  concentration.  Rutgers  [47] 
has  made  an  extensive  survey  and  has  suggested  the  following  correlation  by  Vand  [48, 
49], 

—  =  (1  -  c  - 1 . 1 6c2  )"2  5  (51) 

Mf 

This  correlation  is  derived  for  spherical  particles  and  is  valid  for  both  small  and  large 
particle  concentrations.  In  reality,  the  shape  of  the  encapsulated  particles  might  not  be 
perfectly  spherical  and  hence  the  viscosity  can  be  different.  It  can  be  observed  that  from 
Table- 1  of  [50],  the  measured  experimental  viscosity  values  deviate  from  the  values 
obtained  using  Vand’s  correlation.  This  difference  increases  with  increase  in  particle 
concentration. 


3.4b  Comparison  with  Heat  Transfer  Experiment  Results 

Heat  transfer  experiments  were  performed  in  MMC  heat  sinks  with  an  average  foot  print 
of  1  cm  x  2  cm.  There  are  approximately  441  microchannels  inside  the  heat  sink  and  each 
channel  is  approximately  is  101  pm  wide,  533  pm  high  and  1  mm  long.  MEPCM  slurry 
made  at  BASF  was  used  and  consisted  of  particles  of  size  ranging  from  1  to  5  pm  with  an 
average  diameter  of  4.97  pm  mixed  in  water.  The  core  material  is  n-octadecane  and  the 
shell  material  is  polymethyimetacrylate  (PMMA).  For  numerical  simulation,  viscosity 
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values  provided  in  [50]  were  used.  The  latent  melting  enthalpy  of  PCM  was  measured 
using  differential  scanning  calorimetry  (DSC)  and  was  found  to  be  120  kJ/kg  between  23 
°C  and  29  C.  Properties  used  for  numerical  simulation  are  shown  in  Table  7. 


Table  7:  Properties  of  the  suspension  components  used  for  experiments 


Component 

Density, 

kg/m’ 

Specific 

heat.  J/kg.K 

Thermal  conductivity, 

W/m.K 

Viscosity, 

kg/m.s 

Latent  heat, 

kJ/kg 

Water 

997 

4180 

0.604 

1  x  1  O'3 

- 

MEPCM 

Particle 

867.2 

1899 

0.1643 

- 

120 

Slurry  (10%) 

982.3 

3951 

0.541 

2.3x1 0'-5 

- 

Configuration  1  of  Table  5  is  used  for  numerical  simulation.  Since  the  average  particle 
diameter  is  4.97  pm,  the  same  value  was  used.  Copper  is  used  as  the  wall  material.The 
flow  rate  at  the  microchannel  inlet  was  calculated  based  on  the  number  of  channels  and 
total  flow  rate  at  the  inlet  of  the  heat  sink.  The  wall  thickness  of  the  microchannel  was 
used  such  that  the  base  heat  flux  and  the  total  heat  supplied  is  the  same  as  in  experiments. 
In  other  words,  twice  the  base  area  in  the  simulation  domain  multiplied  by  number  of 
channels  is  equal  to  the  base  are  of  the  heat  sink,  which  is  2  cm  .  Pressure  drop  and  bulk 
temperature  rise  obtained  in  both  cases  are  presented  in  Table  8.  Pressure  drop  predicted 
numerically  ranges  from  8.3  to  13.0%  of  total  pressure  drop  obtained  in  experiments. 
This  suggests  that  the  pressure  drop  inside  manifolds  is  higher  compared  to  that  of 
pressure  drop  inside  the  microchannels,  which  could  be  due  to  poor  manifold  design. 
Proper  design  of  the  heat  sink  can  result  in  more  than  90%  of  pressure  drop  inside 
microchannels  and  less  than  10%  in  manifolds.  Figure  39  shows  the  heat  transfer 
coefficient  obtained  with  numerical  simulation  and  comparison  with  experiment  results. 
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Table  8:  Numerical  results  obtained 


Fluid 

Vhs(GPH) 

Q(W) 

T,n(°C) 

ATbu|k(°C) 

AP  (psi) 

Water 

5 

183.6 

22.1 

8.5 

0.10 

Water 

10 

360.1 

33.8 

8.3 

0.29 

Water 

15 

360.5 

34.2 

5.8 

0.62 

Slurry 

5 

183.6 

22.1 

8.9 

0.18 

Slurry 

10 

360.1 

23.6 

8.7 

0.45 

Slurry 

15 

360.5 

25.2 

6.1 

0.81 

Figure  39:  Comparison  of  numerical  and  experiment  results 


The  predicted  maximum  wall  temperature  using  the  numerical  model  does  not  account 
for  the  thermal  resistance  between  the  heater  surface  and  the  channel  wall.  Hence,  the 
heat  transfer  coefficient  obtained  numerically  is  higher  compared  to  experiment  results. 
Unstructured  mesh  was  created  using  the  default  mesh  option  in  COMSOL  and  the 
results  were  checked  for  grid  independency.  The  results  were  found  to  satisfy  the  global 
energy  balance  at  the  inlet  and  outlet  of  the  channel  w  ithin  4%  balance. 
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Analysis  of  Experiment  Results 


The  performance  of  slurry  is  poorer  compared  to  pure  water  at  all  the  flow  rates.  The 
possible  reasons  for  this  behavior  of  slurry  could  be  due  to: 

•  Low  thermal  conductivity  of  slurry  compared  to  water 

The  resultant  thermal  conductivity  of  slurry  is  lower  compared  to  pure  water  and  hence 
the  heat  transfer  from  wall  to  the  fluid  in  case  of  slurry  is  lower. 

•  Little  or  no  melting  of  PCM  particles  inside  the  channel 

Since  the  length  of  the  channel  is  short  for  the  MMC  channel,  it  is  important  that  the 
PCM  particle  completely  melts  within  its  residence  time.  Figure  40  shows  the  required 
temperature  difference  between  the  fluid  and  the  particle  melting  temperature  for  a  5  pm 
diameter  particle  at  different  velocities.  A  velocity  of  0.25  m/s  corresponds  to  Vhs  =  5 
GPH,  lowest  flow  rate  used  in  the  experiments.  Since  the  flow  is  not  fully  developed 
thermally,  it  is  highly  possible  that  a  large  portion  of  the  fluid  remained  at  the  inlet 
temperature  and  this  required  temperature  difference  could  not  be  achieved  for  all  the 
particles  inside  the  channel. 


Figure  40:  Required  temperature  difference  between  particle  surface  and  PCM 

melting  temperature 
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3.4c  Parametric  Analysis  Results 


Based  on  the  above  results,  it  is  obvious  that  the  performance  of  slurry  depends  on  the 
thermal  conductivity  of  base  fluid.  Also,  the  particle  size  should  be  smaller  so  that  the 
particle  melts  by  the  time  it  exits  the  channel.  Since  most  of  the  fluid  away  from  the  walls 
did  not  absorb  heat  from  the  walls  in  case  of  101  pm  wide  channels,  simulations  were 
also  done  in  channels  of  25  pm  width.  Parametric  study  was  continued  by  varying  the 
base  fluid,  particle  concentration  and  channel  dimensions.  Table  9  shows  the  properties 
of  suspension  components  used. 


Table  9:  Thermophysical  properties  used 


Density 

(kg/m3) 

Specific 

heat 

(J/kg.K) 

Thermal 

conductivity 

(W/m.K) 

Viscosity 

(kg/m.s) 

Latent 

heat 

(J/kg) 

PAO 

783 

2242 

0.143 

4.45x1  O'3 

- 

n-octadecane 

815 

2000 

0.18 

- 

244x10* 

Water 

997 

4180 

0.604 

lxl  0‘3 

- 

Copper 

8700 

385 

400 

- 

- 

Performance  of  water  and  PAO  in  101  pm  wide  channels 

PAO  based  slurry  was  used  for  simulations  in  channels  used  for  experiment,  assuming  a 
particle  diameter  of  100  nm.  The  viscosity  of  the  fluid  was  calculated  using  Equation 
4.26,  since  there  are  no  reported  values  of  dy  namic  viscosity  for  nanoPCM  based  fluids. 
For  comparison,  a  heat  sink  inlet  flow  rate  of  10  GPH  and  a  heat  flux  of  180  W/cm  were 
used  for  both  water  and  PAO  based  fluids.  The  inlet  temperature  used  is  25  JC,  which  is 
within  the  melting  range.  Figure  41  shows  the  bulk  temperature  rise  predicted  for  both 
water  and  PAO.  It  can  be  observed  that  the  bulk  temperature  rise  decreases  with  increase 
in  the  particle  mass  concentration.  Figure  42  shows  the  heat  transfer  coefficient  ratio  (hr  = 
hsiurry/hpuretiuid)  of  both  water  and  PAO  as  base  fluid. 
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Figure  41:  Bulk  temperature  rise  for  water  and  PAO 


Figure  42:  Heat  transfer  coefficient  ratio  for  water  and  PA  O 


In  Figure  42,  hwater  =  61832  W/nr.K  and  hpAo  =  18286  W/m2.K..  The  heat  transfer 
coefficient  decreases  with  increase  in  particle  mass  concentrations  when  water  is  used  as 
the  base  fluid  except  with  mass  concentration  equal  to  0.05.  Slight  increase  in  heat 
transfer  coefficient  at  mass  concentration  of  0.05  indicates  that  the  degradation  in  thermal 
conductivity  is  not  significant  and  hence  the  increase  in  the  specific  heat  compared  to 
base  fluid  helps  the  slurry  performance.  When  PAO  is  used  as  base  fluid,  the  heat  transfer 
coefficient  increased  with  increase  in  concentration.  This  shows  that  for  developing  PCM 
slurry  flows,  thermal  conductivity  plays  a  very  important  role  in  determining  the  slurry 
performance. 
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Performance  in  25  //m  wide  channels 


Since  most  of  the  fluid  did  not  absorb  heat  in  case  of  101  pm  wide  channels  as  the 
thermal  boundary  layer  was  not  developed,  smaller  hydraulic  diameter  was  used  in  order 
to  aid  the  thermal  boundary  to  develop  faster.  Hence  numerical  investigation  was 
continued  by  using  the  simulation  domain  as  configuration  11  (channel  width  25  pm).  The 
total  number  of  channels  inside  the  heat  sink  is  2000.  An  inlet  flow  rate  of  10  GPH  and 
heat  flux  of  180  W/cm2  was  used. 


Figures  43  and  44  show  the  temperature  profile  of  w  ater  based  slurry'  in  both  the  channels 
considered.  It  can  be  observed  that  in  case  of  101  pm  wide  channels,  large  portion  of 
fluid  did  not  absorb  heat  and  in  case  of  25  pm  wide  channels,  the  temperature  profile  is 
more  developed. 


Max:  317.5 


Figure  43:  Temperature  profile  of  water  based  slurry  in  101  pm  wide  channels 
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Figure  44:  Temperature  profile  of  water  based  slurry  in  25  pm  wide  channels 


Figures  45  and  46  show  the  bulk  temperature  rise  and  heat  transfer  coefficient  ratio 
predicted  with  water  as  base  fluid  at  different  fluid  inlet  temperatures.  The  heat  transfer 
coefficient  of  slurry  is  1 .5  times  of  water  at  high  mass  concentration  enabling  a  maximum 
wall  temperature  decrease  of  4  K. 


Figure  45:  Bulk  temperature  rise  in  case  of  25  pm  wide  channels  (base  fluid  - 

water) 
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Figure  46:  Heat  transfer  coefficient  ratio  for  water  (hHttter  =  150695  W/m‘.K) 


Figures  47  and  48  show  the  bulk  temperature  rise  obtained  with  PAO  as  base  fluid  and 
the  heat  transfer  coefficient  of  PAO  slurry  is  tw  ice  of  pure  PAO  at  high  concentration. 


Figure  47:  Bulk  temperature  rise  (base  fluid  -  PAO) 
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Figure  48:  Heat  transfer  coefficient  ratio  ( Base  Fluid  PAO,  lino  =  61930  W/m2.K) 


The  maximum  wall  temperature  for  30%  PAO  slurry  is  15  K  lower  than  the  maximum 
wall  temperature  of  pure  PAO.  Decreasing  the  channel  width  to  25  pm  allow  more 
number  of  fins  within  the  heat  sink.  Compared  to  wider  channels,  more  amount  of  fluid 
absorbs  heat  within  each  channel  and  the  flow'  is  closer  to  thermally  fully  developed  as 
the  thermal  entrance  length  is  short  in  case  of  narrow  channels.  Hence,  using  smaller 
width  channels  enabling  more  heat  transfer  to  the  fluid  helps  to  obtain  better  heat  transfer 
coefficient  for  slurries  compared  to  pure  base  fluid.  Amount  of  PCM  that  participates  in 
heat  absorption  by  the  time  the  fluid  exits  the  channel  depends  on  the  inlet  temperature. 
For  water,  if  the  inlet  temperature  is  around  25.5  °C,  the  heat  transfer  coefficient  is  the 
highest,  where  as  for  PAO,  it  is  23  °C. 
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Ejfecj  ofrnass  concentration  on  pressure  drop 


In  order  to  investigate  the  effect  of  high  concentration  on  pressure  drop,  simulations  were 
run  using  pure  water  and  pure  PAO  to  calculate  Performance  Factor  (PF)  defined  as  in 
Equation  52. 

PF  —  APbase  fluid /APslurry  (52) 


This  parameter  signifies  increase  or  decrease  in  the  pressure  drop  when  a  pure  fluid  is 
used  in  order  to  achieve  the  same  heat  transfer  coefficient  as  slurry.  For  example,  the  heat 
transfer  coefficient  of  water  based  slurry  at  5%  concentration  is  171672  W/irf.K  when 
the  heat  sink  inlet  flow  rate  is  1 0  GPH  and  heat  flux  is  1 80  W/cm  .  In  order  to  achieve  the 
same  heat  transfer  coefficient  with  water,  a  higher  heat  sink  inlet  flow  rate  was  used  and 
the  resultant  pressure  drop  inside  the  channel  was  1.14  times  of  the  pressure  drop  when 
5%  slurry  was  used.  Figure  49  shows  the  PF  with  both  water  and  PAO,  when  the  channel 
width  is  25  pm  and  heat  flux  is  180  W/cm'.  The  heat  sink  inlet  flow  rate  used  for  pure 
water  and  PAO  is  higher  than  10  GPH.  where  as  for  slurry  it  is  10  GPH  for  all  the 
concentrations. 


Figure  49:  Performance  Factor  v.v  cm 

For  high  mass  concentrations,  the  pressure  drop  in  case  of  slurry  is  higher  compared  to 
pure  water  for  same  heat  transfer  coefficient.  It  can  be  concluded  that  particle  mass 
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concentration  of  0.1  has  the  highest  PF  for  both  water  and  PAO  based  slurries.  It  is 
assumed  that  the  microchannel  pressure  drop  dominates  the  total  pressure  drop.  These 
results  are  assumed  true  for  the  parameters  considered  and  might  change  with  flow  rate, 
channel  dimensions,  heat  flux,  fluid  inlet  temperature  etc. 


3.4d  Conclusion 

Thermal  conductivity  of  slurry  plays  a  very  important  role  in  the  cooling  performance  of 
slurry  in  microchannels,  especially  in  manifold  microchannels,  which  provide  flow 
lengths  that  are  comparable  to  the  developing  length  of  the  flow.  Results  show  that  the 
heat  transfer  coefficient  of  water  based  slurry  is  lower  compared  to  pure  water  when  the 
channel  w  idth  is  101  pm.  For  the  same  configuration,  when  the  base  fluid  is  changed  to 
pure  PAO,  that  has  the  thermal  conductivity  equal  to  that  of  the  PCM,  heat  transfer 
coefficient  of  slurry  is  higher  compared  to  pure  PAO  and  the  heat  transfer  coefficient 
increased  with  increase  mass  concentrations.  Thus  in  order  to  achieve  better  performance 
of  slurry  in  developing  flows,  presence  of  PCM  particles  should  enhance  the  thermal 
conductivity  of  the  base  fluid. 


Slurry  performance  depends  on  the  geometric  configuration  of  the  microchannel. 
Numerical  investigation  showed  that  in  microchannels  of  width  25  pm.  water  based 
slurry  too  showed  high  heat  transfer  coefficient  compared  to  pure  water  at  all  mass 
concentrations.  This  is  possibly  because  the  flow  develops  faster  in  smaller  channels. 


When  same  heat  transfer  coefficient  is  desired  with  both  slurry  and  water,  results  showed 
that  higher  mass  concentrations  are  not  favorable  because  of  large  pressure  drop  across 
the  microchannel.  Particle  mass  concentration  of  0.1  showed  the  highest  value  of 
performance  factor  for  the  parameters  considered. 
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Performance  of  PCM  slurries  in  case  of  thermally  developing  flows  depends  on  many 
factors.  Using  a  PCM  that  can  enhance  the  thermal  conductivity  of  slurry,  smaller 
MEPCM  particles  that  can  melt  instantaneously,  narrow  channels  that  can  help  the  flow 
develop  faster,  along  with  right  inlet  temperature  can  help  to  achieve  maximum  benefit  of 
slurry  and  hence  tuning  of  slurry  parameters  is  very  important. 
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SECTION  4:  PLANNED  ACTIVITIES 


The  three  proposed  tasks  are  completed.  In  order  to  analyze  the  slurry  performance  in 
more  detail  and  extend  the  concept  to  a  wider  range  of  parameters,  future  work  may 
involve  the  following: 

•  Synthesis  of  nanoPCM  particles. 

•  Perform  experiments  using  nanoPCM  slurry  in  MMC  heat  sinks. 

•  Find  optimum  configuration  (right  combination  of  NEPCM  particles  and  fluid, 
concentration  depending  on  the  choice  of  suspension  components,  tune  inlet 
temperature  and  velocity)  for  obtaining  the  best  thermal  performance  for  a  given 
application. 
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